$Z$ and $R$ are complex-valued functions with real and imaginary parts in $S_h$ which depend upon $n$, $\ell$ and $j$ and are given by $$\label{E310}
Z=U_{j+1}^{n,1,\ell+1}+{\rm i} \mu U_{j+1}^{n,2,\ell+1},\quad R=r_{j+1}^{n,1,\ell+1}+{\rm i}\mu r_{j+1}^{n,2,\ell+1}.$$

In practice only a finite number of outer and inner iterates are computed at each time step. Specifically, for $i=1,2$, $n\geq 0$, we compute approximations to the outer iterates $U_j^{n,i}$ for $j=1,\ldots, J_{out},$ for some small positive integer $J_{out}$. For each $j$, $0\leq j\leq J_{out}-1$, $U_{j+1}^{n,i}$ is approximated by the last inner iterate $U_{j+1}^{n,i,J_{inn}}$ of the sequence of inner iterates $U_{j+1}^{n,i,\ell}$, $0\leq \ell \leq J_{inn}$ that satisfy linear systems of the form (\[E39\]). $J_{inn}$ and $J_{out}$ are such that $$\left(\sum_{k=1}^2\|U_{j+1}^{n,k,\ell+1}-U_{j+1}^{n,k,\ell}\|^2_{\ell_2}\right)^{1/2}\leq \varepsilon,$$ and $$\left(\sum_{k=1}^2\|U_{j+1}^{n,k}-U_{j}^{n,k}\|^2_{\ell_2}\right)^{1/2}\leq \varepsilon,$$ where $\|v\|_{\ell_2}$ denotes the Euclidean norm of the coefficients of $v\in S_h$ with respect to its basis, and $\varepsilon$ is usually taken to be $10^{-10}$.

Given $U^n$, the required starting values $U_0^{n,i}$ for the outer (Newton) iteration are computed by extrapolation from previous values as $$\label{E311}
U_0^{n,i}=\alpha_{0,i} U^{n}+\alpha_{1,i} U^{n-1}+\alpha_{2,i} U^{n-2}+\alpha_{3,i} U^{n-3},$$ for $i=1,2$, where the coefficients $\alpha_{j,i}$ are such that $U_0^{n,i}$ is the value at $t=t^{n,i}$ of the Lagrange interpolating polynomial of degree at most 3 in $t$ that interpolates to the data $U^{n-j}$ 