even, $k$ odd} \\ X_z(i,j,k,t)=\sqrt{\mu_{i,j,k}}H_z(i,j,k,t), & \mbox{$i$ odd, $j$ odd, $k$
even} \\ Y_x(i,j,k,t)=\sqrt{\varepsilon_{i,j,k}}E_x(i,j,k,t), & \mbox{$i$ odd, $j$ even, $k$ even}
\\ Y_y(i,j,k,t)=\sqrt{\varepsilon_{i,j,k}}E_y(i,j,k,t), & \mbox{$i$ even, $j$ odd, $k$ even} \\
Y_z(i,j,k,t)=\sqrt{\varepsilon_{i,j,k}}E_z(i,j,k,t), & \mbox{$i$ even, $j$ even, $k$ odd}
\\
\end{array} \right.\,,
\label{psidef}$$ for the origin of the coordinate system $(i,j,k)=(0,0,0)$ at the center of the unit cell shown in Fig. \[fig:3D\]. The number of lattice points in the $x$, $y$, and $z$ direction will be denoted by $n_x$, $n_y$, and $n_z$ respectively. As before these numbers are assumed to be odd.

Discretization of the differential operators that appear in Eq. (\[h\]) yields Eq. (\[eqn:mtxeqn\]) in the form $${\frac{\partial}{\partial t}}\Psi(t)={H} \Psi(t)=
\mathop{{\sum}'}_{i=1}^{n_x}
\mathop{{\sum}'}_{j=1}^{n_y}
\mathop{{\sum}'}_{k=1}^{n_z}
\left[
H^{(x)}(i,j,k)+H^{(y)}(i,j,k)+H^{(z)}(i,j,k)
\right]\Psi(t),
\label{3Dsplit}$$ where the superscripts $(x)$, $(y)$ and $(z)$ refer to the derivative with respect to $x$, $y$ and $z$ respectively, $$\begin{aligned}
H^{(x)}(i,j,k)=&+&
\frac{\be^{{\phantom{T}}}_{i,j+1,k}\be^T_{i+1,j+1,k}-
\be^{{\phantom{T}}}_{i+1,j+1,k}\be^T_{i,j+1,k}}
{\delta\sqrt{{\varepsilon}_{i+1,j+1,k}\mu_{i,j+1,k}}}
-
\frac{\be^{{\phantom{T}}}_{i,j,k+1}\be^T_{i+1,j,k+1}-
\be^{{\phantom{T}}}_{i+1,j,k+1}\be^T_{i,j,k+1}}
{\delta\sqrt{{\varepsilon}_{i+1,j,k+1}\mu_{i,j,k+1}}}
\label{hx} \\
&+&
\frac{\be^{{\phantom{T}}}_{i+1,j+1,k}\be^T_{i+2,j+1,k}-
\be^{{\phantom{T}}}_{i+2,j+1,k}\be^T_{i+1,j+1,k}}
{\delta\sqrt{{\varepsilon}_{i+1,j+1,k}\mu_{i+2,j+1,k}}}
-
\frac{\be^{{\phantom{T}}}_{i+1,j,k+1}\be^T_{i+2,j,k+1}-
\be^{{\phantom{T}}}_{i+2,j,k+1}\be^T_{i+1,j,k+1}}
{\delta\sqrt{{\varepsilon}_{i+1,j,k+1}\mu_{i+2,j,k+1}}},
\nonumber\end{aligned}$$ and the expressions for $H^{(y)}(i,j,k)$ and $H^{(z)}(i,j,k)$ follow from Eq. (\[hx\]) by symmetry. Note that we use the triple $(i,j,k)$ to label the $n_xn_yn_z$ unit vectors $\be_{i,j,k}$.

In complete analogy with the 1D and 2D case we split Eq. (\[hx\]) (as well as 