Let us know implement the implicit scheme and show how we can extend the previous algorithm for solving Laplace's or Poisson's equations to the diffusion equation as well. As the reader will notice, this simply implies a slight redefinition of the vector \( \mathbf{b} \) defined in Eq. (20).
To see this, let us first set up the diffusion in two spatial dimensions, with boundary and initial conditions. The \( 2+1 \)-dimensional diffusion equation (with dimensionless variables) reads for a function \( u=u(x,y,t) \) $$ \begin{equation*} \frac{\partial u}{\partial t}= D\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}\right). \end{equation*} $$