 Solution

Choose the centered Dirichlet discrete Laplacian in each direction. With $T=\mathrm{tridiag}(1,-2,1)$ put $A=T\otimes I$, $B=I\otimes T$ and $\mu =\Delta t/h^2$. Thus $h^{-2}A$ and $h^{-2}B$ are the physical derivative approximations; the factor $h^{-2}$ has been absorbed into the Courant number in the written update. Equivalently use the scaled derivative matrices and the time parameter $\Delta t$.

For a one-dimensional vector with endpoint values $v_0=v_{m+1}=0$,Thus $A$ and $B$ are real symmetric negative-definite matrices. An orthonormal eigenvector basis shows that all eigenvalues of each exponential lie in $(0,1]$ for $\mu \ge 0$, soIterating gives $\Vert U^n{\Vert }_2\le \Vert U^0{\Vert }_2$ for every nonnegative step parameter, a mesh-independent unconditional stability bound. Multiplying the squared norm by the grid-area weight $h^2$ leaves the same estimate. In this natural tensor-grid choice $AB=BA=T\otimes T$, so $E(\mu ;A,B)=e^{\mu (A+B)}$: the split update is even the exact semidiscrete diffusion flow. Contractivity of a product of these negative symmetric exponentials would prove stability even without commutation.