[proofplan]
We prove that the discrete error $\|e\|_2 \leq \frac{c}{8} h$ by assembling three estimates. First, we derive the error equation $Ae = \eta$ by subtracting the discrete equation for the exact solution (which includes a truncation error $\eta$) from the discrete equation for the numerical solution. Second, we bound the truncation error vector $\|\eta\|_2 \leq c h^3$ using Taylor's theorem and the grid size relation $m < 1/h$. Third, we bound the inverse $\|A^{-1}\|_2 \leq \frac{1}{8h^2}$ using the explicit eigenvalue formula and the elementary inequality $\sin\theta \geq \frac{2}{\pi}\theta$ for $\theta \in [0, \frac{\pi}{2}]$. Combining these three estimates yields the result.
[/proofplan]
[step:Derive the error equation $Ae = \eta$ from the consistency relation]
Let $\widehat{u}$ denote the exact solution of Poisson's equation on $\Omega = (0,1)^2$, and let $U_{i,j}$ denote the numerical solution at grid point $(x_i, y_j) = (ih, jh)$. Define the error $e_{i,j} := U_{i,j} - \widehat{u}(x_i, y_j)$ for $i, j = 1, \ldots, m$.
The numerical solution satisfies the five-point equation exactly:
\begin{align*}
U_{i-1,j} + U_{i+1,j} + U_{i,j-1} + U_{i,j+1} - 4U_{i,j} = h^2 f_{i,j}.
\end{align*}
By the [Taylor Expansion of the Second Central Difference](/theorems/1363) applied to $\widehat{u} \in C^4(\overline{\Omega})$ in each coordinate direction, the exact solution satisfies
\begin{align*}
\widehat{u}_{i-1,j} + \widehat{u}_{i+1,j} + \widehat{u}_{i,j-1} + \widehat{u}_{i,j+1} - 4\widehat{u}_{i,j} = h^2 \Delta\widehat{u}(x_i, y_j) + \eta_{i,j},
\end{align*}
where $\eta_{i,j}$ is the local truncation error satisfying
\begin{align*}
\eta_{i,j} = \frac{h^4}{12}\left(\frac{\partial^4 \widehat{u}}{\partial x^4}(\xi_i, y_j) + \frac{\partial^4 \widehat{u}}{\partial y^4}(x_i, \zeta_j)\right)
\end{align*}
for intermediate points $\xi_i \in (x_{i-1}, x_{i+1})$ and $\zeta_j \in (y_{j-1}, y_{j+1})$. Since $\Delta\widehat{u} = f$ on $\Omega$, we have $h^2 \Delta\widehat{u}(x_i, y_j) = h^2 f_{i,j}$. Subtracting the exact relation from the numerical one:
\begin{align*}
e_{i-1,j} + e_{i+1,j} + e_{i,j-1} + e_{i,j+1} - 4e_{i,j} = -\eta_{i,j}.
\end{align*}
In matrix form, this reads $Ae = -\eta$. Since $A$ is invertible by the [Symmetry and Negative Definiteness of $A$](/theorems/1367), we obtain $e = -A^{-1}\eta$, and therefore
\begin{align*}
\|e\|_2 \leq \|A^{-1}\|_2 \, \|\eta\|_2.
\end{align*}
[guided]
The idea behind this step is the standard consistency framework: the numerical solution satisfies the discrete equation exactly (by definition), while the exact solution satisfies it only up to a truncation error. Subtracting the two equations isolates the error, which satisfies the same linear system but with the truncation error as the right-hand side.
Let $\widehat{u}$ denote the exact solution of Poisson's equation on $\Omega = (0,1)^2$, and let $U_{i,j}$ denote the numerical solution at grid point $(x_i, y_j) = (ih, jh)$. Define the error $e_{i,j} := U_{i,j} - \widehat{u}(x_i, y_j)$ for $i, j = 1, \ldots, m$.
The numerical solution satisfies the five-point equation exactly:
\begin{align*}
U_{i-1,j} + U_{i+1,j} + U_{i,j-1} + U_{i,j+1} - 4U_{i,j} = h^2 f_{i,j}.
\end{align*}
For the exact solution, we apply the [Taylor Expansion of the Second Central Difference](/theorems/1363) to $\widehat{u} \in C^4(\overline{\Omega})$ in the $x$-direction at the point $(x_i, y_j)$:
\begin{align*}
\widehat{u}(x_i - h, y_j) - 2\widehat{u}(x_i, y_j) + \widehat{u}(x_i + h, y_j) = h^2 \frac{\partial^2 \widehat{u}}{\partial x^2}(x_i, y_j) + \frac{h^4}{12}\frac{\partial^4 \widehat{u}}{\partial x^4}(\xi_i, y_j)
\end{align*}
for some $\xi_i \in (x_{i-1}, x_{i+1})$, and analogously in the $y$-direction. Adding the two one-dimensional expansions:
\begin{align*}
\widehat{u}_{i-1,j} + \widehat{u}_{i+1,j} + \widehat{u}_{i,j-1} + \widehat{u}_{i,j+1} - 4\widehat{u}_{i,j} = h^2 \Delta\widehat{u}(x_i, y_j) + \eta_{i,j},
\end{align*}
where $\eta_{i,j} = \frac{h^4}{12}(\frac{\partial^4 \widehat{u}}{\partial x^4}(\xi_i, y_j) + \frac{\partial^4 \widehat{u}}{\partial y^4}(x_i, \zeta_j))$. Since $\widehat{u}$ solves $\Delta\widehat{u} = f$, the term $h^2 \Delta\widehat{u}(x_i, y_j) = h^2 f_{i,j}$. Subtracting from the numerical equation gives $Ae = -\eta$.
Why the minus sign? The numerical equation has $h^2 f_{i,j}$ on the right and no truncation error, while the exact equation has $h^2 f_{i,j} + \eta_{i,j}$. Subtracting: $(AU)_{i,j} - (A\widehat{u})_{i,j} = h^2 f_{i,j} - (h^2 f_{i,j} + \eta_{i,j}) = -\eta_{i,j}$.
Since $A$ is invertible by the [Symmetry and Negative Definiteness of $A$](/theorems/1367), we obtain $e = -A^{-1}\eta$, and by the submultiplicativity of the operator norm:
\begin{align*}
\|e\|_2 = \|A^{-1}\eta\|_2 \leq \|A^{-1}\|_2 \, \|\eta\|_2.
\end{align*}
The problem has now been reduced to bounding $\|\eta\|_2$ (how well the scheme approximates the PDE locally) and $\|A^{-1}\|_2$ (how much the scheme amplifies errors).
[/guided]
[/step]
[step:Bound the truncation error $\|\eta\|_2 \leq c h^3$]
By the definition of $c$ in the theorem statement, we have
\begin{align*}
|\eta_{i,j}| \leq \frac{h^4}{12}\left(\left|\frac{\partial^4 \widehat{u}}{\partial x^4}(\xi_i, y_j)\right| + \left|\frac{\partial^4 \widehat{u}}{\partial y^4}(x_i, \zeta_j)\right|\right) \leq \frac{h^4}{12} \cdot 12c = c h^4
\end{align*}
for each grid point $(i, j)$. Computing the Euclidean norm of $\eta \in \mathbb{R}^{m^2}$:
\begin{align*}
\|\eta\|_2^2 = \sum_{i=1}^{m} \sum_{j=1}^{m} |\eta_{i,j}|^2 \leq m^2 \cdot c^2 h^8.
\end{align*}
Since $h = \frac{1}{m+1}$, we have $m = \frac{1}{h} - 1 < \frac{1}{h}$, so $m^2 < \frac{1}{h^2}$. Therefore
\begin{align*}
\|\eta\|_2^2 < \frac{1}{h^2} \cdot c^2 h^8 = c^2 h^6,
\end{align*}
and taking square roots:
\begin{align*}
\|\eta\|_2 < c h^3.
\end{align*}
[/step]
[step:Bound $\|A^{-1}\|_2 \leq \frac{1}{8h^2}$ using the explicit eigenvalues]
Since $A$ is real symmetric by the [Symmetry and Negative Definiteness of $A$](/theorems/1367), the spectral theorem gives $\|A\|_2 = \rho(A)$, where $\rho(A) = \max_{k,\ell} |\lambda_{k,\ell}|$ is the spectral radius. The eigenvalues of $A^{-1}$ are $\lambda_{k,\ell}^{-1}$, and $A^{-1}$ is also symmetric, so
\begin{align*}
\|A^{-1}\|_2 = \rho(A^{-1}) = \max_{k,\ell} \frac{1}{|\lambda_{k,\ell}|} = \frac{1}{\min_{k,\ell} |\lambda_{k,\ell}|}.
\end{align*}
By the [Explicit Eigenvalues of $A$](/theorems/1368), $\lambda_{k,\ell} = -4(\sin^2\frac{k\pi h}{2} + \sin^2\frac{\ell\pi h}{2})$, so $|\lambda_{k,\ell}| = 4(\sin^2\frac{k\pi h}{2} + \sin^2\frac{\ell\pi h}{2})$. This is minimised when both $k$ and $\ell$ are as small as possible, i.e., $k = \ell = 1$:
\begin{align*}
\min_{k,\ell} |\lambda_{k,\ell}| = 4 \cdot 2\sin^2\frac{\pi h}{2} = 8\sin^2\frac{\pi h}{2}.
\end{align*}
We apply the elementary inequality $\sin\theta \geq \frac{2}{\pi}\theta$ for $\theta \in [0, \frac{\pi}{2}]$. Since $\frac{\pi h}{2} \in (0, \frac{\pi}{2})$ for $h \in (0, 1)$, we obtain
\begin{align*}
\sin\frac{\pi h}{2} \geq \frac{2}{\pi} \cdot \frac{\pi h}{2} = h.
\end{align*}
Therefore $\sin^2\frac{\pi h}{2} \geq h^2$, and
\begin{align*}
\|A^{-1}\|_2 = \frac{1}{8\sin^2\frac{\pi h}{2}} \leq \frac{1}{8h^2}.
\end{align*}
[guided]
Since $A$ is real symmetric by the [Symmetry and Negative Definiteness of $A$](/theorems/1367), the spectral theorem guarantees that $A$ is unitarily diagonalisable with real eigenvalues. For any real symmetric matrix $B$, the $\ell^2$ operator norm equals the spectral radius: $\|B\|_2 = \rho(B)$. This is because the eigenvectors form an orthonormal basis, and the norm of $Bx$ is maximised by choosing $x$ along the eigenvector with largest eigenvalue magnitude.
Applying this to $A^{-1}$ (which is also symmetric, since the inverse of a symmetric matrix is symmetric):
\begin{align*}
\|A^{-1}\|_2 = \rho(A^{-1}) = \max_{k,\ell} |\lambda_{k,\ell}^{-1}| = \frac{1}{\min_{k,\ell} |\lambda_{k,\ell}|}.
\end{align*}
By the [Explicit Eigenvalues of $A$](/theorems/1368), $|\lambda_{k,\ell}| = 4(\sin^2\frac{k\pi h}{2} + \sin^2\frac{\ell\pi h}{2})$. Since $\sin^2$ is an increasing function on $[0, \frac{\pi}{2}]$, and $\frac{k\pi h}{2} \leq \frac{m\pi h}{2} = \frac{\pi}{2} \cdot \frac{m}{m+1} < \frac{\pi}{2}$ for all valid $k$, the minimum of $|\lambda_{k,\ell}|$ over $k, \ell \in \{1, \ldots, m\}$ is achieved at $k = \ell = 1$:
\begin{align*}
\min_{k,\ell} |\lambda_{k,\ell}| = 8\sin^2\frac{\pi h}{2}.
\end{align*}
To convert this into a power of $h$, we use the concavity inequality $\sin\theta \geq \frac{2}{\pi}\theta$ for $\theta \in [0, \frac{\pi}{2}]$. This inequality holds because $\sin$ is concave on $[0, \frac{\pi}{2}]$ and the line $\frac{2}{\pi}\theta$ connects $\sin(0) = 0$ to $\sin(\frac{\pi}{2}) = 1$; concavity ensures the curve lies above the chord. Applying this with $\theta = \frac{\pi h}{2}$:
\begin{align*}
\sin\frac{\pi h}{2} \geq \frac{2}{\pi} \cdot \frac{\pi h}{2} = h,
\end{align*}
so $8\sin^2\frac{\pi h}{2} \geq 8h^2$, and therefore $\|A^{-1}\|_2 \leq \frac{1}{8h^2}$.
[/guided]
[/step]
[step:Combine the estimates to obtain $\|e\|_2 \leq \frac{c}{8} h$]
Substituting the bounds from the previous two steps into the inequality $\|e\|_2 \leq \|A^{-1}\|_2 \, \|\eta\|_2$:
\begin{align*}
\|e\|_2 \leq \frac{1}{8h^2} \cdot c h^3 = \frac{c}{8} \, h.
\end{align*}
This holds for all sufficiently small $h > 0$ (specifically, for any $h = \frac{1}{m+1}$ with $m \geq 1$), completing the proof that the five-point formula converges with order $O(h)$ in the $\ell^2$ norm.
[/step]