a fully implicit variant of our two-dimensional scheme, with two intermediate arrays, $u^{\star}$ and $u^{\star \star}$,

\[eq:yoshida\] $$\begin{aligned}
u^{\star} & = u(t) + \frac{k}{4} \left( A_{x} u^{\star} + A_{x} u(t) \right)
\label{eq:yoshidaA} \\ 
u^{\star \star} & = u^{\star} + \frac{k}{2} \left( A_{y} u^{\star \star} + A_{y} u^{\star} \right)
\label{eq:yoshidaB} \\
& \qquad \; + \frac{h}{2} \left( (1-u^{\star \star}) u^{\star \star} + (1-u^{\star}) u^{\star} \right)
\nonumber \\ 
u(t+h) & = u^{\star \star} + \frac{k}{4} \left( A_{x} u(t+h) + A_{x} u^{\star \star} \right).
\label{eq:yoshidaC} \end{aligned}$$

In (\[eq:yoshida\]), the operators $A_{x}$ and $A_{y}$ are the same as (\[eq:Amatrix\]) for directions $x$ and $y$, respectively. For simulations on a lattice, $u_{ij}(t) = u(t,x_{0}+(i-1)\Delta
x,y_{0}+(j-1)\Delta y)$, where $1 \leq i \leq N_x,\, 1 \leq j \leq
N_y$ and $\Delta x = \Delta y$, the following gives the action of the $A_{x}, A_{y}$ operators: $$\begin{aligned}
A_{x} u_{i,j} & = & u_{i-1,j} - 2 u_{i,j} + u_{i+1,j}, \\
A_{y} u_{i,j} & = & u_{i,j-1} - 2 u_{i,j} + u_{i,j+1}.\end{aligned}$$ The compression scheme and code outline given in Appendix \[sec:code\] show that only a maximum of one row or column (i.e., max($N_x,N_y)$) of storage is needed for $u^{\star}$ and $u^{\star \star}$.

Again because the fully implicit quadratic vector 