from ours. The Fourier-mode stability analysis performed by Zheng et al.[@Zheng1] does not generalize to the case of spatially varying permittivity and permeability whereas in our approach, the algorithms are unconditionally stable by construction.

Our presentation is organized as follows: The basic theoretical concepts are given in Sec. \[sec2\]. In Sec. \[sec3\] we explain the general philosophy that underlies the Suzuki approach to construct algorithms that are unconditionally stable. In Sec. \[sec4\], we show in detail how to implement these ideas for the case of the TDME in one, two and three spatial dimensions, using algorithms that are accurate up to second and fourth order in the time step. In Sec. \[sec5\] we explain how we analyze the data generated by these algorithms. In Sec. \[sec6\] we present the results of numerical simulations for the physical systems that we selected as examples to test the algorithms. Our conclusions are given in Sec. \[sec7\].

Theory {#sec2}
======

The model system we consider in this paper describes EM fields in a $d$-dimensional ($d=1,2,3$) medium with spatially varying permittivity and/or permeability, surrounded by a perfectly conducting box. In the absence of free charges and currents, the EM fields in such a system satisfy 