Understanding the matrix form of Finite Difference Scheme I am trying to implement a numerical finite difference central difference method to solve the elliptic equation 
$u_{xx} + u_{yy} = sin(\pi y)(2-6x-(\pi x)^2(1-x) )$ for a 3x2 grid. 
This yields the equations: 
$u_{xx} = (u_{i+1, j} + u_{i-1,j} - 2u_{i,j}) + \bigg( \frac{hx^2}{hy^2} \bigg) (u_{i,j+1} + u_{j-1} - 2u_{i,j})   = h^2_{x} sin(\pi y_{j}) (2-6x_{i} -(\pi x_{i})^2 (1-x_{i}) )$
I am trying to understand why the linear system $A u = b$ is constructed. 
The exercise gives u as $u = (u11,u12,u21,u22,u31,u32)^T$ which I can understand as its a 3x2 grid. 
But I do not understand why A looks like this (note $\alpha = \frac{h_{x}}{h_{y}}$) 
$ A = \begin{pmatrix}
-2(1+\alpha) & \alpha & 1 & 0 & 0 & 0 \\
\alpha & -2(1+\alpha) & 0 & 1 & 0 & 0 \\
1 & 0 & -2(1+\alpha) & 1 & \alpha & 0 \\
0 & 1 & \alpha & -2(1+\alpha) & 0 & 1 \\
0 & 0 & 1 & 0 & -2(1+\alpha) & \alpha\\
0 & 0 & 0 & \alpha & 1 & -2(1+\alpha)
\end{pmatrix}
$
Note that $b = h^2_{x} sin(\pi y_{j}) (2-6x_{i} -(\pi x_{i})^2 (1-x_{i}) )$. 
I have tried to substitute for i = 1 and j = 1 for $u_{11}$ on the lhs but that gave some coefficient of u01 which does not make sense. 
How can I arrive at each of the 36 coefficients of $A_{ij} $ from the data I am given ? This is probably very easy but I cannot see it :/  
Edit: Additional information: 
Spatial domain taken to be a unit square $(x,y) \in [0,1]^2 $ with $h_{x} = \frac{1}{4}$ and $h_{y} = \frac{1}{3}$. And we have a spatial grid with 3 interior points and 2 interior points in y. Boundary conditions will be imposed as $ u_{0,j} = 0, u_{4,j} = 0 , u_{i,0} = 0, u_{i,3} = 0 $
Edit: Additional FD to Matrix form examples 
Its not just this matrix I am also having troublinf arriving at A for $(2.6)$ page 4, $2.8$ page 9, $2.12 $ page 10, etc ... on the following notes https://www.uni-muenster.de/imperia/md/content/physik_tp/lectures/ws2016-2017/num_methods_i/advection.pdf. Any hints/tips for a calculation are welcomed! 
 A: We have that 
$$
u_{xx} + u_{yy} \approx \frac{u_{i+1,j} - 2u_{i,j} + u_{i-1,j}}{h_x^2} + \frac{u_{i,j+1} - 2u_{i,j} + u_{i,j-1}}{h_y^2} = \sin(\pi y_j)(2-6x_i-(\pi x_i)^2(1-x_i)).
$$
Multiplying by $h_x^2$ and defining $\alpha = h_x/h_y$ gives
$$
u_{i+1,j} - 2u_{i,j} + u_{i-1,j}+ \alpha^2(u_{i,j+1} - 2u_{i,j} + u_{i,j-1})= \underbrace{h_x^2\sin(\pi y_j)(2-6x_i-(\pi x_i)^2(1-x_i))}_{=b_{i,j}}.
$$
Define the right-hand side to be $b_{i,j}$. Then, rearranging, we get
\begin{align}\tag{1}
u_{i+1,j} + u_{i-1,j} + (-2-2\alpha^2)u_{i,j}+\alpha^2u_{i,j+1} + \alpha^2u_{i,j-1} = b_{i,j}.
\end{align}
Suppose that we stack all the unknowns $u_{i,j}$ into a big column vector
$$
\mathbf{u} = \begin{bmatrix} u_{1,1} \\ u_{1,2} \\ u_{1,N_y} \\ u_{2,1} \\ \vdots \\ u_{N_x,N_y} \end{bmatrix},
$$
where there are $N_x$ grid points in the $x$ direction and $N_y$ in the $y$ direction. Let's write Equation (1) in matrix form, this shall correspond to just one row of the matrix $A$. If we then stack all these rows together we shall obtain the matrix $A$. Observe that Equation (1) is equivalent to
\begin{equation}\tag{2}
\begin{bmatrix} \cdots & 1 & \cdots & \alpha^2 & -2 - 2\alpha^2& \alpha^2 & \cdots & 1 & \cdots \end{bmatrix} \underbrace{\begin{bmatrix} \vdots \\ u_{i-1,j} \\ \vdots \\ u_{i,j-1} \\ u_{i,j} \\ u_{i,j+1} \\ \vdots \\ u_{i+1,j} \\ \vdots \end{bmatrix}}_{\mathbf{u}} = b_{i,j}.
\end{equation}
All the entries in the row not explicitly written out are zero. Note that the "$1$
"'s is are total of $N_y$ to the left and right of the $-2-2\alpha^2$ entry. (Understanding Check 1: Can you explain why?)
There is one issue with what we wrote out: what if our node is on the boundary (if $i = 1$ or $N_x$ or $j = 1$ or $N_y$)? The solution is simple in this case simply don't include the entry in the row vector if it is on the boundary. This corresponds to enforcing zero on the boundary of the domain. (Understanding Check 2: Can you explain why?)
Equation (2) corresponds to enforcing the discretized PDE at the $(i,j)$th position. Thus, if we stack the rows from Equation (2) in the order $(1,1)$, $(1,2)$, ..., $(N_x,N_y)$, same as we stacked $\mathbf{u}$, we get the matrix $A$. In your case with $N_x = 3$ and $N_y = 2$, we get 
\begin{equation}
A = \begin{bmatrix} -2-2\alpha^2 & \alpha^2 & 1 & 0 & 0 & 0 \\
\alpha^2 & -2-2\alpha^2 & 0 & 1 & 0 & 0 \\
1 & 0 & -2-2\alpha^2 & \alpha^2 & 1 & 0 \\
0 & 1 & \alpha^2 & -2-2\alpha^2 & 0 & 1 \\
0 & 0 & 1 & 0 & -2-2\alpha^2 & \alpha^2 \\
0 & 0 & 0 & 1 & \alpha^2 & -2-2\alpha^2\end{bmatrix}
\end{equation}
and $A\mathbf{u} = \mathbf{b}$. (Understanding Check 3: How should $\mathbf{b}$ be constructed to make this work?)
As you can see, the matrix $A$ has some very distinctive patterns. (Understanding Check 4: See if you can write $A$ if instead $N_x = 2$ and $N_y= 3$?)

Response to First Comment The biggest conceptual challenge when writing down the matrix of a finite difference scheme is that the matrix is two-dimensional object, indexed by a row $I$ and a column $J$. However, each index $I$ and $J$ represents a location of a grid point $(i_x,i_y)$ or $(j_x,j_y)$ in the 2D domain on which your solving the underlying PDE. So in a sense, the matrix is more like a four-dimensional object, where $a_{I,J}$ represents the coefficient of $u_{j_x,j_y}$ in Equation (2) enforced at the gridpoint $u_{i_x,i_y}$. More simply, each row $I$ of the equation corresponds to enforcing the PDE at a gridpoint $(i_x,i_y)$.
The question is how to get between $I$ and $i_x$ and $i_y$. Since we traverse the grid in the order $(1,1)$, $(1,2)$, ..., $(1,N_y)$, $(2,1)$,..., $(N_x,N_y)$, with some algebra we can convince ourselves that the point $(i_x,i_y)$ corresponds to row $I = N_y(i_x-1) + i_y$ of the matrix. This gives us a way of converting the point $(i_x,i_y)$ in the PDE discretization to the row (or column) $I$ in the matrix. (Conversely, $i_x$ and $i_y$ can be obtained by dividing $I$ by $N_y$ and computing the quotient and the remainder. The details are somewhat yucky but they can be worked out.)
To work out $a_{I,J}$ we first determine which gridpoints $(i_x,i_y)$ and $(j_x,j_y)$ correspond to the row $I$ and column $J$. We then have three cases:


*

*If $i_x = j_x$ and $i_y = j_y$, then $a_{I,J} = -2 - 2\alpha^2$.

*If $i_x = j_x$ and $i_y = j_y-1$ or $i_y = j_y+1$, then $a_{I,J} = \alpha^2$.

*If $i_y = j_y$ and $i_x = j_x-1$ or $i_x = j_x+1$, then $a_{I,J} = 1$.


(Understanding Check 5: How do I get my three rules by looking at Equation (1)?)

Response to Second Comment Your question is about how to get $i_x$, $i_y$, $j_x$, and $j_y$ from $I$ and $J$. It's important to note that there's an unwritten fourth case I didn't include in my last update.


*If none of the other three cases apply, $a_{I,J} = 0$.


Now how do you get $(i_x,i_y)$ from $I$? There's an easy way and a hard way. The easy way is to count. In your example we have $N_x = 3$ and $N_y = 2$ so $\mathbf{u} = (u_{11},u_{12},u_{21},u_{22},u_{31},u_{32})$. $I=2$ corresponds to $(i_x,i_y) =(1,2)$ since $u_{12}$ is the second element of $\mathbf{u}$. $J=3$ corresponds to $(j_x,j_y) = (2,1)$. Therefore, since we're not in any of the three other cases in my answer, we're in case four and thus $a_{I,J} = a_{23} = 0$.
"Just count" doesn't seem very mathematical. You might want an actual formula for $i_x$ and $i_y$ given $I$. These correspond to the "yucky details" I mentioned in my last post. Let me flesh them out now.
Let $\text{mod}(m,n)$ and $\text{div}(m,n)$ represent the remainder and quotient when the integer $n$ is divided into $m$. For example $\text{mod}(11,3) = 2$ and $\text{div}(11,3) = 3$. We see that $m = n\times \text{div}(m,n) + \text{mod}(m,n)$. Thus, since we know $I = N_y(i_x-1) + i_y$, we have $I - 1 = N_y(i_x-1) + (i_y-1)$, so we have $i_x - 1 = \text{div}(I-1,N_y)$ and $i_y -1 = \text{mod}(I-1,N_y)$. Solving, we have
\begin{align}
i_x &= \text{div}(I-1,N_y)+1, \\
i_y &= \text{mod}(I-1,N_y)+1.
\end{align}
