the second-order central-difference approximation to the first derivative with respect to $x$, we obtain $$\begin{aligned}
{\frac{\partial}{\partial t}}X_y(i,t) & = &
\frac{1}{\delta\sqrt{\mu_i}} \left(
\frac{Y_z(i+1,t)}{\sqrt{\varepsilon_{i+1}}}-
\frac{Y_z(i-1,t)}{\sqrt{\varepsilon_{i-1}}}\right),
\label{eqn:discrTMx} \\
{\frac{\partial}{\partial t}}Y_z(j,t) & = &
\frac{1}{\delta\sqrt{\varepsilon_j}} \left(
\frac{X_y(j+1,t)}{\sqrt{\mu_{j+1}}}-
\frac{X_y(j-1,t)}{\sqrt{\mu_{j-1}}}\right),
\label{eqn:discrTMy}\end{aligned}$$ where the integer $i$ labels the grid points and $\delta$ denotes the distance between two next-nearest neighbor lattice points (hence the absence of a factor two in the nominator). For notational simplicity we will, from now on, specify the spatial coordinates through the lattice index $i$, e.g. $X_y(i,t)$ stands for $X_y(x=(i+1)\delta/2,t)$.

(8.0,1.6) (0.5,0)

Following Yee [@Yee66] it is convenient to assign $X_y(i,t)$ and $Y_z(j,t)$ to the odd, respectively, even numbered lattice site, as shown in Fig. \[fig:fig1\] for a grid of $n$ points. The equations (\[eqn:discrTMx\]) and (\[eqn:discrTMy\]) can now be combined into one equation of the form Eq. (\[TDME\]) by introducing the $n$-dimensional vector $$\Psi(i,t) =
\left\{ \begin{array}{lll} X_y(i,t)=\sqrt{\mu_i}H_y(i,t), & \mbox{$i$ odd}
\\ Y_z(i,t)=\sqrt{\varepsilon_i}E_y(i,t), & \mbox{$i$ even}
\end{array} \right..
\label{eqn:eqnindx}$$ The vector $\Psi(t)$ describes both the magnetic and the electric field on the lattice points $i=1,\ldots,n$. As usual the $i$-th element of $\Psi(t)$ is given by the inner product $\Psi(i,t)=\be^T_i\cdot\Psi(t)$ where $\be_i$ denotes the $i$-th unit vector in the $n$-dimensional vector space. Using 