initial state (at $t=0$) of the EM fields satifies Eqs.(\[max3\]) and (\[max4\]), time-integration of the TDME using $U_k(\tau)$ yields a solution that will not satisfy Eqs.(\[max3\]) and (\[max4\]). However, for algorithm $U_k(\tau)$ the deviations from zero vanish as $\tau^k$ so that in practice these errors are under control and can be made sufficiently small for practical purposes.

Three dimensions {#subsec43}
----------------

In terms of $\Psi(t)=(\bX(t),\bY(t))^T$ for a 3D system, Eq. (\[eqn:mtxeqn\]) reads $${\frac{\partial}{\partial t}} \Psi(t) = {\mathcal H} \Psi(t) =
\left(\begin{array}{cc} {\phantom{+}}0^{\phantom{T}} &
{\phantom{+}}{h}^{\phantom{T}} \\
-h^T & {\phantom{+}}0^{\phantom{T}} \end{array}\right)\Psi(t),\label{eqn:mtx3D}$$ where $h$ is given by $$h=\left(\begin{array}{ccc}
0 & {\phantom{+}}\frac{1}{{\sqrt{\mu}}}{\frac{\partial}{\partial z}}\frac{1}{{\sqrt{\varepsilon}}}
& -\frac{1}{{\sqrt{\mu}}}{\frac{\partial}{\partial y}}\frac{1}{{\sqrt{\varepsilon}}} \\
-\frac{1}{{\sqrt{\mu}}}{\frac{\partial}{\partial z}}\frac{1}{{\sqrt{\varepsilon}}} & 0 &
{\phantom{+}}\frac{1}{{\sqrt{\mu}}}{\frac{\partial}{\partial x}}\frac{1}{{\sqrt{\varepsilon}}} \\
{\phantom{+}}\frac{1}{{\sqrt{\mu}}}{\frac{\partial}{\partial y}}\frac{1}{{\sqrt{\varepsilon}}} &
-\frac{1}{{\sqrt{\mu}}}{\frac{\partial}{\partial x}}\frac{1}{{\sqrt{\varepsilon}}} & 0
\end{array}\right).
\label{h}$$ We discretize the spatial coordinates by adopting the standard Yee grid [@Yee66]. We show a unit cell of this grid in Fig. \[fig:3D\].

(7,6.2) (0.5,0)

In analogy to the systems in 1D and 2D, we assign the EM fields to the lattice points such that the boundary conditions Eq. (\[boundcond\]) are automatically satisfied. The elements of the vector $\Psi(i,j,k,t)$ are given by $$\Psi(i,j,k,t) = \left\{ \begin{array}{lll} X_x(i,j,k,t)=\sqrt{\mu_{i,j,k}}H_x(i,j,k,t), &
\mbox{$i$ even, $j$ odd, $k$ odd}
\\ X_y(i,j,k,t)=\sqrt{\mu_{i,j,k}}H_y(i,j,k,t), & \mbox{$i$ odd, $j$
