<s>[INST] 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)
$$
 [/INST] 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.
</s>