$H^{(y)}(i,j,k)$ and $H^{(z)}(i,j,k)$) in two parts and obtain for the first-order approximation to $U(\tau)$ $$\begin{aligned}
U_1(\tau) & = &
e^{\tau{H_{1}^{(x)}}}
e^{\tau{H_{2}^{(x)}}}
e^{\tau{H_{1}^{(y)}}}
e^{\tau{H_{2}^{(y)}}}
e^{\tau{H_{1}^{(z)}}}
e^{\tau{H_{2}^{(z)}}}\,,
\label{U1D3}\end{aligned}$$ where for instance $$H_{1}^{(z)} =
\mathop{{\sum}'}_{i=1}^{n_x}
\mathop{{\sum}'}_{j=1}^{n_y}
\mathop{{\sum}'}_{k=1}^{n_z}
\left(
\frac{\be^{{\phantom{T}}}_{i,j+1,k}\be^T_{i,j+1,k+1}-
\be^{{\phantom{T}}}_{i,j+1,k+1}\be^T_{i,j+1,k}}
{\delta\sqrt{{\varepsilon}_{i,j+1,k+1}\mu_{i,j+1,k}}}-
\frac{\be^{{\phantom{T}}}_{i+1,j,k}\be^T_{i+1,j,k+1}-
\be^{{\phantom{T}}}_{i+1,j,k+1}\be^T_{i+1,j,k}}
{\delta\sqrt{{\varepsilon}_{i+1,j,k+1}\mu_{i+1,j,k}}}
\right).
\label{h1z}$$ Note that each contribution to Eq. (\[h1z\]) acts on a different pair of elements of $\Psi(t)$. Hence each of the matrix exponentials in Eq. (\[U1D3\]) acts on one quarter of all the lattice points. Performing the time-step operation Eq. (\[U1D3\]) involves only two sweeps over all lattice points.

By construction the algorithm defined by Eq. (\[U1D3\]) is unconditionally stable and so are the higher-order algorithms defined by $U_2(\tau$) and $U_4(\tau)$. Each contribution to e.g. Eq. (\[h1z\]) is of the form Eq. (\[h1\]) and hence its matrix exponential can be calculated in exactly the same manner as in the 1D case (see Eq. (\[twobytwo\])). The divergence of the EM fields in 3D is, for the same reason as in the 2D case, not conserved but decreases as $\tau^k$.

Implementation: Summary {#subsec44}
-----------------------

The notation required to write down the algorithms in mathematical form might give the impression that these algorithms are difficult to program. Actually that is not 