Solution (source code)

= Solution

For the Fourier mode $u_m^n=g^ne^{im\xi}$, substitution gives
$$
\left(1-\frac{i\mu}2\sin\xi\right)g=1+\frac{i\mu}2\sin\xi,\qquad \boxed{g=\frac{1+i\mu\sin\xi/2}{1-i\mu\sin\xi/2}.}
$$
The denominator is never zero for real $\mu,\xi$, and numerator and denominator are complex conjugates, so $|g|=1$. Parseval's identity therefore makes each step norm-preserving on square-summable grid data. \b[The method is von Neumann stable for every $\mu>0$], with no upper Courant restriction. It is an implicit nondissipative scheme; unconditional <stability> does not remove phase error or guarantee accuracy at a large time step.