Solution (source code)

= Solution

The $[1/1]$ <Padé approximant> of the exponential is $r(z)=(1+z/2)/(1-z/2)$. Replacing each exponential in its original order gives the split <Crank-Nicolson method>
$$
u^{n+1}=C_A(\delta)C_B(\delta)u^n,
\qquad C_A(\delta)=(I-\delta A/2)^{-1}(I+\delta A/2),
$$
with the analogous definition for $B$. Operationally, first solve the $B$ substep and then the $A$ substep. The matrices $I-\delta A/2$, $I-\delta B/2$ are invertible for every $\delta\ge0$, since $A,B$ are symmetric negative definite.

An orthonormal eigenbasis of $A$ diagonalizes $C_A$. Its <eigenvalue> for $\lambda\le0$ is
$$
r(\delta\lambda)=\frac{1+\delta\lambda/2}{1-\delta\lambda/2},\qquad
|r(\delta\lambda)|\le1.
$$
Therefore $\|C_A\|_h\le1$, and similarly $\|C_B\|_h\le1$. Submultiplicativity gives the <contractivity of split Crank-Nicolson diffusion>:
$$
\boxed{\|u^{n+1}\|_h\le\|u^n\|_h\quad\text{for every }\mu=\Delta t/h^2\ge0.}
$$
No commutativity or simultaneous diagonalization is needed for this product bound. This proves unconditional discrete $L^2$ stability, not positivity of every stencil weight or unconditional maximum-norm monotonicity. Nor do the individually second-order rational substeps make the ordered splitting second order: when $[A,B]\ne0$, its leading splitting defect remains $\delta^2[A,B]/2$.