Solution (source code)

= Solution

Split the generator as $iA-iV$. One <Strang splitting> step is
$$
\boxed{
\mathbf u^{n+1}
=e^{-ihV/2}e^{ihA}e^{-ihV/2}\mathbf u^n
}.
$$
The two potential half-steps are componentwise multiplications. The circulant matrix $A$ is diagonalized by the <discrete Fourier transform>; its eigenvalues are
$$
\lambda_j=-4M^2\sin^2\left(\frac{\pi j}{2M}\right),
\qquad j=0,\ldots,2M-1.
$$
Thus the middle step consists of a <Fast Fourier transform>, multiplication of Fourier coefficient $j$ by $e^{ih\lambda_j}$, and an inverse transform. Its cost is $O(M\log M)$ and every factor is unitary, so the implementation preserves the discrete norm exactly up to roundoff.