Confusion with numerical method A book that I'm studying says this: 
For a particular value $n$ consider the matrix $A$ of dimension $n^2\times n^2$ with all its elements zero except: $$a_{i,i}=4+h^2$$
$$a_{i,i+1}=a_{i,i-1}=a_{i,i+n}=a_{i,i-n}=-1$$
with $h=1/(n+1)$.
The system $Au=b$ where $b_i=4h^2$ for every $i$ offers an approximation of the solution $u(x,y)$ of the PDE $$\frac{\partial ^2u}{\partial x^2} + \frac{\partial ^2u}{\partial y^2}-u=4$$ with $u(x,y)=0$ in the region $0\leq x,y\leq 1$.
I don't understand because $u$ is an $n^2$ dimension vector and $u(x,y)$ is  a function. Any help on what this approximation gives?
 A: Let $\Omega \subset \mathbb{R}^2$ denote the unit square, i.e., $$\Omega  = (0,1) \times (0,1). $$ Our overall objective is to solve the partial differential equation
$$ \Delta u(x,y) + \alpha u(x,y) = f(x,y), \quad (x,y) \in \Omega$$
subject to homogeneous boundary conditions, i.e.,
$$ u(x,y) = 0, \quad (x,y) \in \partial \Omega.$$
The original problem concerns the case of $$\alpha = -1, \quad f(x,y) = 4.$$
We shall use a space central approximation of the 2nd derivative twice. Specifically, if $v \in C^2(\overline{\Omega},\mathbb{R})$, then
$$ \frac{\partial^2 v}{\partial x^2}(x,y) = \frac{v(x+h,y)-2 v(x,y) + v(x-h,y)}{h^2} + O(h^2), \quad h \rightarrow 0, \quad h \not = 0$$
and
$$ \frac{\partial^2 v}{\partial y^2}(x,y) = \frac{v(x,y+h)-2 v(x,y) + v(x,y-h)}{h^2} + O(h^2), \quad h \rightarrow 0, \quad h \not = 0.$$
These statement are true for all $(x,y) \in \Omega$. We now move to cover the closed unit square $$\bar{\Omega} = [0,1] \times [0,1]$$ with a uniform grid. We shall use the same spacing $h$ between the grid point in either direction. We choose a positive integer $m$ and select $h$ such that $mh = 1$. Our grid consist of the points $$(x_i,y_j) = (ih,jh), \quad (i,j) \in \{0,1,2,\dots,m\} \times \{0,1,2,\dotsc,m\}.$$
We mention in passing that there are $(m-1)^2$ grid points that are contained in the open set $\Omega$, specifically the grid points
$$ (x_i,y_j), \quad (i,j) \in \{1,2,\dotsc,n\} \times \{1,2,\dotsc,n\}, \quad n = m - 1.$$
This completes our preparations and we are now ready to tackle the main issues, namely how to store and compute approximations of the solution $u$ at the internal grid points. Since $u$ is zero on the boundary of $\Omega$ there is no compelling reason for storing this information and we are free to concentrate on the $n^2$ internal grid points. Suppose that we are given $n^2$ numbers $v_{ij}$ such that $$v_{ij} \approx u(x_i,y_j), \quad (i,j) \in \{1,2,\dotsc,n\} \times \{1,2,\dotsc,n\}$$
is a good approximation. How do we store this data? Certainly, there is more than one choice, but we choose to define a matrix $V \in \mathbb{R}^{n \times n}$ by $V = [v_{ij}]$. This choice is appealing if you imagine that you are standing at the point $(m,m)$ and are looking towards $(0,0)$. You will find that the axes align nicely with the rows and columns of the matrix. We are now at the first critical point of the analysis. We define a tridiagonal matrix $A \in \mathbb{R}^{n \times n}$ as follows, $$ h^2 a_{ij} = \begin{cases} -2 & i=j \\ 1 & |i-j| = 1 \\ 0 & \text{otherwise}\end{cases}.$$ This is the celebrated discrete Laplacian. It is straightforward, but important to verify that
$$(AV)_{ij} = \frac{v_{i-1,j} - 2v_{ij} + v_{i+1,j}}{h^2}.$$
We also have
$$(VA)_{ij} = \frac{v_{i,j-1}- 2v_{ij} + v_{i,j+1}}{h^2}.$$
These statements are true at all internal grid points provided that we set $v_{ij} = 0$ on the boundary nodes.
Now let $u \in C^2(\bar{\Omega},\mathbb{R})$ denote the solution of our differential equation with homogeneous boundary conditions. Let $U = [u_{ij}] \in \mathbb{R}^{n \times n}$ denote the matrix given by $$u_{ij} = u(x_i,y_j).$$
Then by virtue of the aforementioned expansions we have
$$ (AU + UA + \alpha U)_{ij} - f_{ij} =  O(h^2), \quad h \rightarrow 0, \quad h \not = 0,$$
where $F = [f_{ij}] \in \mathbb{R}^{n \times n}$ is given by $$f_{ij} = f(x_i,y_j).$$
This is the reason why we try to solve the matrix equation
$$ AV + VA + \alpha V = F$$
with respect to $V = [v_{ij}] \in \mathbb{R}^{n \times n}$ and use the components of $V$ to approximate $u$, hoping that $$v_{ij} \approx u(x_i,y_j)$$ is a good approximation when $h$ is sufficiently small.
We are now at the second critical point of the analysis. We recognize that the map $\phi : \mathbb{R}^{n \times n} \rightarrow \mathbb{R}^{n \times n}$ given by $$\phi(V) = AV + VA + \alpha V$$ is linear in the input variable $V$. Now let $\text{vec} : \mathbb{R}^{n \times n} \rightarrow \mathbb{R}^{n^2}$ be the linear map such that $\text{vec}(V)$ is formed by stacking the columns of $V$ on top of each other. This is not precise mathematical definition, but it reflects the column major format that is normally used to store dense matrices. If we wish to be fanatical about it, then we use the equation $$ (\text{vec}(V))_{(j-1)n+i} = v_{ij}$$
to define $\text{vec}(V)$. With this definition, we have
$$ AV + VA + \alpha V = F$$
if and only if
$$ (I \otimes A + A \otimes I + \alpha I \otimes I) \text{vec}(V) = \text{vec}(F)$$
Here $I$ is the $n$-by-$n$ identity matrix and $\otimes$ denotes the Kronecker product, see this link for the basic definition and properties. The nontrivial but elementary property that is critical in our context is the identity $$\text{vec}(AXB) = (B^T \otimes A) \text{vec}(X).$$
Our new equation is a standard linear system with a coefficient matrix that is $n^2$-by-$n^2$ and a right-hand side that is a vector with $n^2$ components. The unknown vector has $n^2$ components. There are two differences between the OP's system and the system that has been described here. One is trivial scaling by $h^2$, while the other is due to a sign error. The OP's system is consistent with the differential equation $$- \Delta u + u = 4$$ but this is a minor issue that is unrelated to the conceptual problems that where created when the instructor overloaded the notation without ample explanation and used $u$ to denote both the smooth solution and the discrete approximation.
