Finite difference method applied to the 2D time-independent Schrödinger equation The 1D Schrödinger equation
$$
E\psi=\left[-\frac{\hbar^2}{2m}\frac{d^2}{dx^2}+\hat{V}\right]\psi
$$
can be solved using a finite difference scheme. The kinetic energy operator is represented by the matrix
$$
-\frac{\hbar^2}{2m}\left( \begin{array}{cccc}
2 & -1 & 0  & \dots & \dots & 0  \\
-1 & 2 & -1 & \dots & \dots & 0 \\
0 & -1 & 2 & -1 & \dots & 0 \\
\vdots & \vdots & \vdots & \ddots &  & \vdots \\
\vdots & \vdots & \vdots &  & \ddots & -1 \\
0 & \dots & \dots & 0 & -1 & 2 \\
\end{array} \right)
$$
Similarly, the potential energy operator has the form
$$
\left( \begin{array}{cccc}
V_1 & 0 & 0  & \dots & \dots & 0  \\
0 & V_2 & 0 & \dots & \dots & 0 \\
0 & 0 & V_3 & 0 & \dots & 0 \\
\vdots & \vdots & \vdots & \ddots &  & \vdots \\
\vdots & \vdots & \vdots &  & \ddots & 0 \\
0 & \dots & \dots & 0 & 0 & V_N \\
\end{array} \right)
$$
My question is: what is the general form for the 2D case?
I know the above matrices are generated from the three-point stencil, where
$$
\frac{\partial^2 \psi}{\partial x^2} = \frac{\psi_{i+1}+\psi_{i-1}-2\psi_{i}}{\Delta x^2}
$$
So with a five-point stencil for the 2d case, I'll have
$$
\frac{\partial^2\psi}{\partial x^2} + \frac{\partial^2\psi}{\partial y^2}= \frac{\psi_{i+1,j}+\psi_{i-1,j}-2\psi_{i,j}}{\Delta x^2} + \frac{\psi_{i,j+1}+\psi_{i,j-1}-2\psi_{i,j}}{\Delta y^2}
$$
But I'm unsure of how to convert this into the relevant matrices.
[edit] I think I've managed to derive a solution by following this video.
The kinetic energy operator is
$$
-\frac{\hbar^2}{2m}\left( \begin{array}{cccc}
4 & -1 & 0  & -1 & 0 & 0 & 0 & 0 \\
-1 & 4 & -1 & 0 & -1 & 0 & 0 & 0 \\
0 & -1 & 4 & 0 & 0 & -1& 0 & 0\\
-1 & 0 & 0 & 4 & -1 & 0& -1 & 0\\
0 & -1 & 0 & -1 & 4 & -1 & 0 & 0\\
0 & 0 & -1 & 0 & -1 & 4 &  0 & 0\\
 &  &  & \vdots &  & &  & \\
\end{array} \right)
$$
And the potential energy operator is
$$
\left( \begin{array}{cccc}
v_{11} & 0 & 0  & 0 & 0 & 0 & 0 & 0 \\
0 & v_{21} & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & \ddots & 0 & 0 & 0& 0 & 0\\
0 & 0 & 0 & v_{N1} & 0 & 0& 0 & 0\\
0 & 0 & 0 & 0 & v_{22} & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 0 & \ddots &  0 & 0\\
 &  &  & \vdots &  & &  & v_{NM} \\
\end{array} \right)
$$
 A: You'll want to use Kronecker sum of discrete Laplacians for 2D or higher-dimensional Laplacian. If your potential energy operator is separable into sum of 1D functions, then you can similarly use Kronecker sum for potential energy operator. Otherwise you just have to appropriately sample your potential energy into corresponding matrix so as to match the layout of your discrete Laplacian.
A: The natural solution to solve this problem is to store $\psi$ into a 2D matrix $\Psi = \Psi_{ij}$.
Let's call
$$D = \begin{pmatrix}
  2 & -1 & 0 & \cdots & \cdots & 0 \\
  -1 & 2 & -1 & \cdots & \cdots & 0 \\ 
  0 & -1 & 2 & -1 & \cdots & 0 \\
  \vdots & \vdots & \vdots  & \ddots & & \vdots\\
  \vdots & \vdots & \vdots & & \ddots & -1 \\
  0 & 0 & 0 & \cdots & -1 & 2 
\end{pmatrix}$$
then, as you noticed, $D\Psi$ represents $\frac{\textrm{d}^2\psi}{\textrm{d}x^2}$.
To find a way of calculating $\frac{\textrm{d}^2\psi}{\textrm{d}y^2}$, you can think about the fact that $\Psi^\top$ represents $\psi$, where the $x$ and $y$ axis are swapped, so $D\Psi^\top$ represents $\frac{\textrm{d}\psi}{\textrm{d}y}$, with swapped axis. However, $D$ is symetric, so it's judicious to write $D=\, D^\top$, which finally gives $D\Psi^\top = (\Psi D)^\top$: we can now swap the axis once again, so $\frac{\textrm{d}^2\psi}{\textrm{d}x^2}$ is represented by $D\Psi$.
Then, you can represent $\Delta\psi$ by $D\Psi + \Psi D$
Notice that here I implicitly assumed that $\Psi$ was a square matrix: if it's not the case, you will have to create to matrix $D_x$ and $D_y$, which are similar to $D$, but with different dimensions.
