Implicit 2D finite difference linear system I am familiar with the 1D implicit method which solves the heat equation with homogeneous Dirichlet conditions,
$$u_i^{k-1} = \big(1+2\lambda)u_i^k - \lambda \big( u_{i+1}^k +  u_{i-1}^k \big) $$  which can be rearranged in to a linear system $$\underline{u} ^{k-1}  = A \underline{u} ^{k}$$ where $$A = \begin{pmatrix}
 1+2\lambda & -\lambda & 0 & \dots & 0\\ 
 -\lambda & 1+2\lambda & -\lambda & \ddots & \vdots\\ 
 0 & -\lambda & \ddots & \ddots & 0\\
 \vdots & \ddots & \ddots & \ddots & -\lambda\\
0 & \dots & 0 & -\lambda & 1+2\lambda\\
 \end{pmatrix}.$$
However, I am becoming increasingly confused about what the linear system for the 2 Dimensional case would look like. I am aware that the so called algorithm would now be,
$$  u_{i,j}^{k-1} = \big(1+4\lambda)u_{i,j}^k - \lambda \big( u_{i+1,j}^k +  u_{i-1,j}^k + u_{i,j+1}^k +  u_{i,j-1}^k \big). $$
$\textbf{My question - What is the linear system for this algorithm?}$
I know that it can have the same format as the 1D case but I'm unsure what the vector $\underline{u}$ would be since there are now two different spatial dimensions. (i.e we have both $i$ and $j$ varying). I am also unsure of what the coefficient matrix would be. Thanks in advanced for any help.
 A: Basically you have just to make one decision: How do you want to order the vector $\boldsymbol{u}$. That can be done in row-major ordering, i.e.,
$$\boldsymbol{u} = \begin{pmatrix} u_{1, 1}& u_{2, 1} & \dots u_{N_x, 1} & u_{1, 2} & u_{2, 2} & \dots & u_{N_x, N_y} \end{pmatrix}^T $$
or in column-major
$$\boldsymbol{u} = \begin{pmatrix} u_{1, 1}& u_{1, 2} & \dots u_{1, N_y} & u_{2, 1} & u_{2, 1} & \dots & u_{N_x, N_y} \end{pmatrix}^T $$
The qualitative structure of the matrix is then the same. You have two more diagonals (one below and one above the main diagonal). The "distance" from the diagonal element to the "far" off diagonals is then determined by your ordering.
For row-major, the distance in $\boldsymbol{u}$ from $u_{i,j} $ to $u_{i, j \pm 1}$ is $\pm N_x$, while for column-major, the distance from $u_{i,j} $ to $u_{i \pm 1, j }$ is $\pm N_y $.
Thus, the matrix A reads:
$$A = \begin{pmatrix}
 1+2\lambda & -\lambda & \underbrace{0}_{\{N_x, N_y\} - 2} & -\lambda & 0 & \dots & 0 \\ 
 -\lambda & 1+2\lambda & -\lambda & \underbrace{0}_{\{N_x, N_y\} - 2} -\lambda & 0 & \dots & 0\\ 
 0 & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots &\\
-\lambda & \underbrace{0}_{\{N_x, N_y\} - 2} & -\lambda  &
 1+2\lambda & -\lambda \underbrace{0}_{\{N_x, N_y\} - 2} & -\lambda & 0 \\
0 &  & \dots & -\lambda & \underbrace{0}_{\{N_x, N_y\} - 2} & -\lambda  &
 1+2\lambda
 \end{pmatrix}.$$
