In two dimensions, our vector of eigenstates becomes the array
$$ \begin{bmatrix} \psi_{11} & \psi_{12} & \dots & \psi_{1N}\\ \psi_{21} & \psi_{22} & \dots & \psi_{2N}\\ \vdots\\ \psi_{N1} & \psi_{N2} & \dots & \psi_{NN} \end{bmatrix} $$
This can be represented in vector form by
$$ \begin{bmatrix} \psi_{11}\\ \psi_{12}\\ \vdots\\ \psi_{1N}\\ \psi_{21}\\ \psi_{22}\\ \vdots\\ \psi_{2N}\\ \vdots\\ \psi_{N1}\\ \psi_{N2}\\ \vdots\\ \psi_{NN} \end{bmatrix} $$
which has a length of $N^2$. Consider for simplicity the case where $N=3$. The form of $\frac{\partial^2}{\partial{x^2}}$ is
$$ \frac{1}{\Delta x^2} \begin{bmatrix} -2 & 1 & 0\\ 1 & -2 & 1\\ 0 & 1 & -2 \end{bmatrix}\equiv D. $$
In two dimensions, the array of eigenstates is condensed into a vector as shown previously. We define the second derivative using the Kronecker Product:
$$ I\otimes D=\frac{1}{\Delta{x}^2} \begin{bmatrix} D & 0 & 0\\ 0 & D & 0\\ 0 & 0 & D \end{bmatrix} $$
you can check that this takes the second derivative of the condensed vector form.
Similarly, the second derivative in the $y$-direction is
$$ D\otimes I = \frac{1}{\Delta{y^2}} \begin{bmatrix} -2I & I & 0\\ I & -2I & I\\ 0 & I & -2I \end{bmatrix} $$
Therefore,
$$ \frac{\partial^2}{\partial{x^2}}+\frac{\partial^2}{\partial{y^2}}=I\otimes D+D\otimes I=D\oplus D $$
where $\oplus$ is the Kronecker Sum. The two dimensional Schrodinger Equation becomes
$$ \bigg[-\frac 12(D\oplus D)+m\Delta{x^2}VI\bigg]\psi=m\Delta{x^2}E\psi $$
where $V$ is the matrix containing the potentials.
Plotting the probability distribution for $|\psi_{0}|^2$: