Solution (source code)

= Solution
{c}

Applying the <implicit midpoint rule> to the <simple harmonic oscillator> system gives
$$
y_{n+1}-y_n=\frac h2(z_{n+1}+z_n),\qquad
z_{n+1}-z_n=-\frac h2(y_{n+1}+y_n).
$$
Factor the change in the <quadratic form> $y^2+z^2$ and substitute these equations:
$$
\begin{aligned}
&(y_{n+1}^2+z_{n+1}^2)-(y_n^2+z_n^2)\\
&=(y_{n+1}-y_n)(y_{n+1}+y_n)
+(z_{n+1}-z_n)(z_{n+1}+z_n)\\
&=\frac h2(z_{n+1}+z_n)(y_{n+1}+y_n)
-\frac h2(y_{n+1}+y_n)(z_{n+1}+z_n)=0.
\end{aligned}
$$
Thus every step preserves the <first integral>
$$
\boxed{y_{n+1}^2+z_{n+1}^2=y_n^2+z_n^2
=y_0^2+z_0^2.}
$$
For completeness, the implicit equations are solvable for every real $h$: putting $s=h/2$, their coefficient <matrix> has <determinant> $1+s^2>0$ and the update is
$$
\begin{pmatrix}y_{n+1}\\z_{n+1}\end{pmatrix}
=\frac1{1+s^2}
\begin{pmatrix}1-s^2&2s\\-2s&1-s^2\end{pmatrix}
\begin{pmatrix}y_n\\z_n\end{pmatrix}.
$$
This is an <orthogonal matrix>, as follows from $(1-s^2)^2+4s^2=(1+s^2)^2$. It therefore preserves the <Euclidean norm>, explaining geometrically the exact energy conservation. This is a concrete instance of <Runge-Kutta conservation of quadratic invariants>. Preservation of the circle does not assert that the numerical and exact solutions have the same phase at a fixed nonzero step size.