What is the discrete (finite difference) Laplacian operator $\Delta_h$ in $\mathbb{R}^N$? What is the corresponding matrix $A_h$ such that $\Delta_h f = A_h \cdot f$?
On Wikipedia and on most textbooks only one or two-dimensional case is considered: $$\Delta_h f(x,y)\approx {\frac {f(x-h,y)+f(x+h,y)+f(x,y-h)+f(x,y+h)-4f(x,y)}{h^{2}}}.$$
On
There are many different possible finite difference approximations of many different orders of accuracy. In numerical methods a common choice is to use ''centered in space:'' $$f(x+h) = f(x) + f'(x)h + f''(x)\frac{h^2}{2} + f'''(x)\frac{h^3}{6} + f''''(x)\frac{h^4}{24}$$
$$f(x-h) = f(x) - f'(x)h + f''(x)\frac{h^2}{2}- f'''(x)\frac{h^3}{6}+f''''(x)\frac{h^4}{24}$$
$$\Longrightarrow\frac{f(x+h) + f(x-h) - 2f(x)}{h^2}= f''(x) + f''''(\xi)\frac{h^2}{12} $$ $$\Longrightarrow f''(x) = \frac{f(x+h) + f(x-h) - 2f(x)}{h^2} + \mathcal{O}(h^2)$$
This will correspond to a differentiation matrix: $$D_{+-x}f(x) = \frac{1}{h^2}\begin{pmatrix} &-2 & 1& &\\ &1 &\ddots &\ddots & \\ & &\ddots & \ddots & 1 \\ & & & 1 & -2 \end{pmatrix}\begin{pmatrix}f_1(x) \\ f_2(x) \\ \vdots \\ f_{n-1}(x) \\ f_k(x) \end{pmatrix}$$ where $k$ is total number of grid points used
what you are looking at is the analogous idea for a 2 variable function
In general, the formula for $f(x_1,\dots,x_n)$ will be:
$$\nabla^2 f = \sum_{j=1}^n \frac{\partial^2 f}{dx_j^2} = \left(\sum_{j=1}^n D_{+-x_j}\right)f +\mathcal{O}(h^2)$$ (assuming you use the same grid spacing for each variable)
In one dimension we have the approximation for the second derivative:
$$f^{(2)}(x)\approx \frac{f(x+h) + f(x-h) -2 f(x)}{h^2}$$
This then gives you the approximation for the second partial derivative in some particular direction, the Laplacian is the sum of the second derivatives w.r.t. all the coordinates.