$$
\newcommand{\uex}{{u_{\small\mbox{e}}}}
\newcommand{\half}{\frac{1}{2}}
\newcommand{\tp}{\thinspace .}
\newcommand{\Oof}[1]{\mathcal{O}(#1)}
\newcommand{\x}{\boldsymbol{x}}
\newcommand{\dfc}{\alpha} % diffusion coefficient
\newcommand{\Ix}{\mathcal{I}_x}
\newcommand{\Iy}{\mathcal{I}_y}
\newcommand{\If}{\mathcal{I}_s} % for FEM
\newcommand{\Ifd}{{I_d}} % for FEM
\newcommand{\basphi}{\varphi}
\newcommand{\baspsi}{\psi}
\newcommand{\refphi}{\tilde\basphi}
\newcommand{\xno}[1]{x_{#1}}
\newcommand{\dX}{\, \mathrm{d}X}
\newcommand{\dx}{\, \mathrm{d}x}
\newcommand{\ds}{\, \mathrm{d}s}
$$
Newton's method: the Jacobian and its sparsity
$$ J_{i,j,r,s} = \frac{\partial F_{i,j}}{\partial u_{r,s}} $$
Newton system:
$$ \sum_{r\in\Ix}\sum_{s\in\Iy}J_{i,j,r,s}\delta u_{r,s} = -F_{i,j},
\quad i\in\Ix,\ j\in\Iy\tp$$
But \( F_{i,j} \) contains only \( u_{i\pm 1,j} \),
\( u_{i,j\pm 1} \), and \( u_{i,j} \). We get nonzero contributions
only for
\( J_{i,j,i-1,j} \), \( J_{i,j,i+1,j} \), \( J_{i,j,i,j-1} \), \( J_{i,j,i,j+1} \),
and \( J_{i,j,i,j} \). The Newton system collapses to
$$
\begin{align*}
J_{i,j,r,s}\delta u_{r,s} =
J_{i,j,i,j}\delta u_{i,j} & +
J_{i,j,i-1,j}\delta u_{i-1,j} +\\
& J_{i,j,i+1,j}\delta u_{i+1,j} +
J_{i,j,i,j-1}\delta u_{i,j-1}
+ J_{i,j,i,j+1}\delta u_{i,j+1}
\end{align*}
$$