Jacobi's algorithm extended to the diffusion equation in two dimensions
For future time steps, only the boundary values are determined
and we need to solve the equations for the interior part in an iterative way similar to what was done for Laplace's or Poisson's equations.
To see this, we rewrite the previous equation as
$$
\begin{equation*}
u^{l}_{i,j}= \frac{1}{1+4\alpha}\left[\alpha(u^{l}_{i+1,j}+u^{l}_{i-1,j}+u^{l}_{i,j+1}+u^{l}_{i,j-1})+u^{l-1}_{i,j}\right],
\end{equation*}
$$
or in a more compact form as
$$
\begin{equation}
\tag{22}
u^{l}_{i,j}= \frac{1}{1+4\alpha}\left[\alpha\Delta^l_{ij}+u^{l-1}_{i,j}\right],
\end{equation}
$$
with \( \Delta^l_{ij}= \left[u^l_{i,j+1}+u^l_{i,j-1}+u^l_{i+1,j}+u^l_{i-1,j}\right] \).
This equation has essentially the same structure as Eq.
(19), except that
the function \( \rho_{ij} \) is replaced by the solution at a previous time step \( l-1 \). Furthermore, the diagonal matrix elements
are now given by \( 1+4\alpha \), while the non-zero non-diagonal matrix elements equal \( \alpha \). This matrix is also positive definite, meaning in turn that
iterative schemes like the Jacobi or the Gauss-Seidel methods will converge to the desired solution after a given number of iterations.