Solution (source code)

= Solution

The zero <Dirichlet boundary conditions> permit an exact finite-interval calculation. Let $T=\operatorname{tridiag}(1,-2,1)$ on the $M$ interior grid points. The <Dirichlet discrete Laplacian> has an orthogonal <eigenvector> basis
$$
v_m^{(j)}=\sin\frac{jm\pi}{M+1},\qquad
Tv^{(j)}=-4s_jv^{(j)},\qquad
s_j=\sin^2\frac{j\pi}{2(M+1)},\quad 1\leq j\leq M.
$$
Consequently the update matrix $G=(I-rT)^{-1}$ has modal amplification factors
$$
g_j=\frac1{1+4rs_j}.
$$
For every physical $r\geq0$, all denominators are at least one. <Orthogonal diagonalization> therefore proves
$$
\|G^nU^0\|_d\leq\|U^0\|_d,\qquad
\|U\|_d^2=d\sum_{m=1}^M|U_m|^2,
$$
uniformly in $M$, $r$ and $n$. The same bound controls initial perturbations; a forcing increment is propagated by contraction, so successive increments accumulate at most by their sum. \b[The highest-order consistent method is unconditionally stable for all $r\geq0$.] This proof incorporates the boundary conditions, whereas a periodic <Fourier mode> calculation alone would not do so.

For completeness, if “range” is interpreted algebraically to include negative $r$ on a fixed grid, the exact power-stable range is
$$
\boxed{r\in[0,\infty)\ \cup\ \left(-\infty,-\frac1{2\sin^2(\pi/[2(M+1)])}\right].}
$$
Indeed, for $r<0$ every denominator is less than one, so $|g_j|\leq1$ requires $1+4rs_j\leq-1$ for every $j$. The strongest restriction comes from $s_1$. Equality gives a simple amplification factor $-1$ and is allowed because the update is orthogonally diagonalizable. Other negative values either amplify a mode or make the update singular. This extra branch describes backward time stepping on a fixed spatial grid; it is not a positive-time diffusion discretization, and no fixed negative $r$ remains in that branch as $M\to\infty$.