Jacobi's algorithm extended to the diffusion equation in two dimensions

We assume that we have a square lattice of length \( L \) with equally many mesh points in the \( x \) and \( y \) directions. Setting the diffusion constant \( D=1 \) and using the shorthand notation \( u_{xx}={\partial^2 u}/{\partial x^2} \) etc for the second derivatives and \( u_t={\partial u}/{\partial t} \) for the time derivative, we have, with a given set of boundary and initial conditions, $$ \begin{equation*} \begin{array}{cc}u_t= u_{xx}+u_{yy}& x, y\in(0,L), t>0 \\ u(x,y,0) = g(x,y)& x, y\in (0,L) \\ u(0,y,t)=u(L,y,t)=u(x,0,t)=u(x,L,t)0 & t > 0\\ \end{array} \end{equation*} $$