Solution (source code)

= Solution

Use the corrected same-node method identified above, with $\delta=1/M$, interior indices $1,\ldots,M-1$ and zero endpoint values. Its spatial <matrix> $L_\delta$ is real symmetric. Summation by parts gives, for a vector extended by zero at the endpoints,
$$
-v^TL_\delta v=\delta^{-2}\sum_{j=0}^{M-1}a_{j+1/2}(v_{j+1}-v_j)^2.
$$
Hence $-L_\delta$ is positive definite, and $\sum_j(v_{j+1}-v_j)^2\leq4\sum_jv_j^2$ bounds its <eigenvalues> by $4a_+/\delta^2$. The <Forward Euler method> has amplification <matrix> $I+kL_\delta$. Its <eigenvalues> lie in $[1-4\mu a_+,1]$, so their moduli are at most one when
$$
\boxed{0<\mu\leq\frac1{2a_+}.}
$$
Symmetry makes this a contraction in the mesh-weighted <discrete L2 norm>, uniformly over all grids and all coefficient samples satisfying the bounds. There is also a <maximum norm> proof: the update weights are $\mu a_{m-1/2}$, $1-\mu(a_{m-1/2}+a_{m+1/2})$, and $\mu a_{m+1/2}$. They are nonnegative and sum to one, giving a contraction after the boundary values are included.

This mesh-independent bound is sharp. For $a\equiv a_+$, the spatial eigenvectors are $\sin(j\pi m/M)$ and their amplification factors are
$$
1-4\mu a_+\sin^2\frac{j\pi}{2M},\qquad j=1,\ldots,M-1.
$$
If $\mu>1/(2a_+)$, sufficiently fine grids have a highest-frequency factor below $-1$, giving exponential growth in the step number.

On one fixed grid, rather than uniformly under refinement, the sharp universal <stability of a numerical method> bound in the <discrete L2 norm> is slightly larger:
$$
\mu\leq\frac1{2a_+\cos^2(\pi/(2M))}.
$$
Indeed the quadratic form above is bounded by $a_+$ times the constant-coefficient difference form, whose largest <eigenvalue> is $4\cos^2(\pi/(2M))/\delta^2$; equality is realized by the constant coefficient. Its limit as $M\to\infty$ is the boxed bound. No such fixed-grid <matrix> conclusion repairs the ill-defined shifted-boundary update in the literal printed formula.