Solution (source code)

= Solution

Use <von Neumann stability analysis> on the infinite grid. For the <Fourier mode> $u_{m,k}=\widehat u(t)e^{i(m\xi+k\eta)}$, the semidiscrete <eigenvalue> is
$$
\lambda_h(\xi,\eta)=\frac1h\left[-\frac32-\frac12e^{-i\xi}+\frac12e^{i\xi}
+2e^{i\eta}-\frac12e^{2i\eta}\right]
=\frac{-(1-\cos\eta)^2+i[\sin\xi+(2-\cos\eta)\sin\eta]}h.
$$
Its real part is nonpositive for every frequency. The mode evolves by $\widehat u(t)=e^{t\lambda_h}\widehat u(0)$, whose modulus does not grow. The <discrete Fourier transform> and <Parseval identity> therefore yield
$$
\boxed{\|u(t)\|_{\ell_h^2}\le\|u(0)\|_{\ell_h^2}\quad(t\ge0),}
$$
with a bound independent of $h$, where $\|u\|_{\ell_h^2}^2=h^2\sum_{m,k}|u_{m,k}|^2$. Hence \b[the semidiscrete method is stable]. Frequencies with $\eta=0$ are nondissipative rather than unstable.

The same sign can be seen directly by <summation by parts>: the centred $x$ operator is skew-adjoint, and for the unitary $y$ shift $S_y$,
$$
\operatorname{Re}\langle u,D_{y,+,2}u\rangle_h
=-\frac1{4h}\|(S_y-I)^2u\|_h^2\le0.
$$
This is an energy version of the <dissipative second-order forward advection stencil>. A subsequent time integrator would need its own stability condition; semidiscrete stability does not assert that every fully discrete method is stable.