Solution (source code)

= Solution

Write $s=\mu_x+\mu_y$. The update can be rearranged as
$$
\begin{aligned}
(1+s)U_{ij}^{m+1}-\frac{\mu_x}2(U_{i+1,j}^{m+1}+U_{i-1,j}^{m+1})-\frac{\mu_y}2(U_{i,j+1}^{m+1}+U_{i,j-1}^{m+1})
={}&(1-s)U_{ij}^m+\frac{\mu_x}2(U_{i+1,j}^m+U_{i-1,j}^m)+\frac{\mu_y}2(U_{i,j+1}^m+U_{i,j-1}^m).
\end{aligned}
$$
When $s\leq1$, the right side is a <convex combination> of old values: its coefficients are nonnegative and sum to one. Assume all old values lie in $[U_{\min},U_{\max}]$. A new maximum $v>U_{\max}$ is interior because boundary values are zero and $U_{\max}\geq0$. At its node the left side is at least $(1+s)v-sv=v$, while the right side is at most $U_{\max}$, a contradiction. A new minimum $v<U_{\min}$ similarly makes the left side at most $v$ and the right side at least $U_{\min}$, again a contradiction.

The initial values satisfy the bounds by construction; <mathematical induction> proves
$$
\boxed{U_{\min}\leq U_{ij}^m\leq U_{\max}\quad\text{at every node and time level}.}
$$
The new-time matrix is a <strictly diagonally dominant matrix> with positive diagonal and nonpositive off-diagonal entries, so the linear update is uniquely solvable. This <Crank-Nicolson maximum principle under a step restriction> concerns preservation of value bounds, which is a stronger requirement than scalar <A-stability>.