matrices and the vectors of length two, and overwrite the same two elements.

It also follows immediately that performing a time step with algorithms based on Eq. (\[U1D\]) takes ${\cal O}(K)$ plane rotations where $K$ is the total number of elements of the vector $\Psi(t)$ (which is less or equal to the number of grid points). This renders the algorithm efficient: The number of operations to complete one time step scales linearly with the number of grid points. Also note that there is a high degree of intrinsic parallelism in this class of algorithms. In principle, the $(n-1)/2$ matrix-vector multiplications that implement $e^{\tau H_{1}}$ or $e^{\tau H_{2}}$ can be done in parallel.

In the absence of external currents (see below), updating the EM field values of a 3D system using the the Yee algorithm requires 6 arithmetic operations (see Eq.(33) in Ref.[@Taflove]) whereas the second-order algorithm $U_2(\tau)$ requires 33 arithmetic operations. For 1D and 2D problems the ratio is 9/4 and 21/6 respectively. Thus, in terms of CPU time, the price paid for the unconditional stability of the algorithms is not that much and for some applications (see below) may well be worth paying.

An important aspect that we have 