Understanding the matrix equation used to solve 2-D Poisson Equation with non-uniform grid I am working on a module in my class that is meant to teach us about Jacobi/Gauss-Seidel methods of solving matrix equations in conjunction with Poisson's equation, and am having trouble with setting up my matrix.
$$\nabla^2 \phi = 4 \pi \rho$$
with finite difference stencil:
$$\frac{\phi_{i-1,j} - 2\phi_{i,j} + \phi_{i+1,j}}{\Delta x^2} + \frac{\phi_{i,j+1} -2\phi_{i,j} + \phi_{i,j+1}}{\Delta y^2} = 4\pi\rho_{i,j}$$
The goal is to transform this stencil into a matrix equation in two different ways. Once with Dirichlet boundary conditions and once with periodic boundary conditions.
$$ 
\phi_{0,j} = \phi_{N,j}, \forall j\\
\phi_{i, 0} = \phi_{i,N}, \forall i\\
$$
and for Direchlet on the boundaries:
$$
\phi_{i,j} = 1 
$$
I successfully programmed an algorithm using both Jacobi or Gauss Seidel method to solve a matrix, and now I am trying to understand how to put this finite difference stencil into a matrix. I decided to tackle the case of periodic boundary conditions first:
My work so far
Since I am specifically told not to assume $\Delta x = \Delta y$, I rearranged the above equation to find my matrix coefficients with the following substitutions: 
$$
a = 1/\Delta x^2 \\
b = 1/\Delta y^2 \\
c = a + b\\
$$
gives
$$
\frac{1}{4\pi}[b\phi_{i,j-1} + a\phi_{i-1,j} - 2c\phi_{i,j} + a\phi_{i+1,j} + b\phi_{i,j+1} ]
$$
My best attempt at arranging this equation into a matrix looks like:
$$
\mathbf{A} =
\begin{bmatrix}
    -2c & a & b & 0 & \dots &  b & a \\
    a & -2c & a & b & 0 & \dots & b \\
    b & a & -2c & a & b & 0 & \dots  \\
    0 & \ddots & \ddots & \ddots & \ddots &   &  \vdots \\
     \vdots & \\
    b & 0 & \dots & b & a & -2c & a\\
    a & b & 0 & \dots & b & a & -2c \\
\end{bmatrix} \\
$$
Where the solution vector must be ordered as follows:
$$
\begin{align}
& (i_{min}, j_{min})\\
& (i_{min} + 1, j_{min}) \\
& (i_{min} + 2, j_{min}) \\
& \vdots \\
& (i_{max}, j_{min})\\ 
& (i_{min}, j_{min} +1) \\
& (i_{min} + 1, j_{min} + 1) \\
& \vdots \\
& (i_{max} - 1, j_{max}) \\
& (i_{max}, j_{max})
\end{align}
$$
Note: This is my first matrix in TeX, so forgive me if it is a bit ugly. 
So, the matrix equation $\frac{ab}{4\pi}\mathbf{A}\mathbf{\phi} = \mathbf{\rho}$ should yield a system of equations that makes sense to the finite difference stencil, right? 
Well when I test it out, I get about half of what I'm looking for and then some stuff that doesn't make sense even when considering periodic boundary conditions. For example I will expand the matrix equation for $\rho_{2,0}$:
$$
\rho_{2,0} = \frac{1}{4\pi}[b \phi_{(0,0)} + a\phi_{(1,0)} -2c\phi_{(2,0)} + a\phi_{(3,0)} + b\phi_{(4,0)}]
$$
but this should, by the finite difference stencil be:
$$
\rho_{2,0} = \frac{1}{4\pi}[b \phi_{(2,-1)} + a\phi_{(1,0)} -2c\phi_{(2,0)} + a\phi_{(3,0)} + b\phi_{(2,1)}]
$$
So it looks almost right, but no matter how I rearrange my matrix some terms are wrong. I have tried multiple arrangements, but this one gets the most correct terms when I start expanding. 
Where am I going wrong?
 A: Some remarks that are to large for a comment, but not really an answer. 
An important information is missing, but I assume that you are working on the unit square $Ω=I_x×I_y=[0,1]^2$. And $I_x$ is divided in equidistant subintervals $$0=t^x_0≤t^x_1≤…≤t^x_{m^x+1}=1,$$
with $t^x_{i+1}-t^x_{i}=Δx$, and $I_y$ is divided in the same way.  
(I personally like the notation with $m$ representing the inner points and $0$ and $m+1$ are boundary points, but use whatever you like best.)
So your grid looks something like this: 
(0,1)                (1,1)
   15--16--17--18--19 
    |   |   |   |   |    
   10--11--12--13--14
    |   |   |   |   |
    5---6---7---8---9
    |   |   |   |   |
    0---1---2---3---4 
(0,0)                (1,0)

The numbering you introduced is often called "lexicographical numbering" in the literature. Here in this example $m_x = 3$, $m_y=2$. The matrix size is $n=(m_x+2)(m_y+2)$. 
Extracting $\frac{ab}{\cdot}$ from the equation does not give you anything, except a lot of confusion. You wrote this: 
$$\frac{ab}{4\pi}[a^{-1}\phi_{i,j-1} + b^{-1}\phi_{i-1,j} - 2\phi_{i,j}(a^{-1} + b^{-1}) + b^{-1}\phi_{i+1,j} + a^{-1}\phi_{i,j+1}] $$
which is nothing else than
$$\frac{1}{4\pi}[b\phi_{i,j-1} + a\phi_{i-1,j} - 2\phi_{i,j}(b + a) + a\phi_{i+1,j} + b\phi_{i,j+1}]  $$
which looks more like the original stencil, as $b$ is simply the distance between two points in $y$-direction. And that is indicated by a change in the $j$-index.
I don't want to introduce a new notation here, but if I were you I would have used $Δ\tilde{x}=\frac{1}{Δx^2}$ to save space (if necessary).
About the stencil: The stencil represents the matrix row of a node in the centre of the grid. For nodes at the boundary  you have the boundary conditions (see below.)
E.g. for row $k:=7$, representing point 7 you have  
(0,1)                (1,1)
   15--16--17--18--19 
    |   |   |   |   |    
   10--11--12--13--14
    |   |   ∥   |   |
    5---6===7===8---9
    |   |   ∥   |   |
    0---1---2---3---4 
(0,0)                (1,0)

which results in: 
$$\frac{1}{4\pi}[b\phi_{2} + a\phi_{6} - 2\phi_{7}(b + a) + a\phi_{8} + b\phi_{12}]  $$
This will give you the matrix row: 
$$\begin{matrix} … & 0 &b & 0&…&0& a & 2(a+b)&a & 0 & … & 0 &b & 0&… \\ 
  &  &2 & &&& 6 & 7&8 &  &  &  &12 & &  \\
  &  &k-(m_x+2) & &&& k-1 & k&k+1 &  &  &  &k+(m_x+2) & & \end{matrix}$$
First row: Value
Second row: column-index for the specific row 7
Third row: column-index for an arbitrary row $k$
Zero-Boundary-Condition
As I stated above the matrix-size is $n=(m_x+2)(m_y+2)$, but for zero boundary conditions you can reduce the system to $\tilde{n}=m_xm_y$.
E.g. the equation for point 2 reads: 
$$ φ_{3,0}=0 $$
Thus this translates to a row: 
$$\left[\begin{matrix} 0& …& 0 & 1 & 0 & … &0\end{matrix}\right]=0$$
This row does not add any value to the matrix system and the row can be erased from the system. And as the zero boundary condition also introduces a $0$ in the row for node 7 above: 
$$\left[\begin{matrix} … & 0 &\color{red}{0} & 0&…&0& a & 2(a+b)&a & 0 & … & 0 &b & 0&… \end{matrix}\right]$$
You can also erase column 7 from the system. 
I hope that this is can be understood.
Periodic Boundary Condition
You did not state how exactly your boundary condition looks like, but this could be something like:
$$φ(x,0) = φ(x,1) $$
which can be formed to 
$$φ(x,0) - φ(x,1) = 0,$$
which again can be represented in the equation system.
