Development of solution of  Del^2 U(x,y) = C(x,y) 
boundary value PDE by solving simultaneous equations.

From the second order numerical difference 
d^2u/dx^2 + d^2u/dy^2 = 
   1/h^2(U(-1h,0)+U(1h,0)+U(0,-1h)+U(0,1h)-4U(0,0)) +O(h^2) = C(x,y)
or U(-1h,0) + U(1h,0) + U(0,-1h) + U(0,1h) -4U(0,0) = C(x,y)*h^2


Using subscripts for a matrix layout
solve for i=1,nx-1  j=1,ny-1
 U(i-1,j) + U(i+1,j) + U(i,j-1) + U(i,j+1) -4U(i,j) = C(i,j)*h^2

The solution U is not known, except at the boundary, B(0,j), B(nx,j)
B(i,0) and B(i,ny)

A matrix layout for a 5 by 5 matrix is

     +--------+--------+--------+--------+--------+
     | U(4,0) | U(4,1) | U(4,2) | U(4,3) | U(4,4) |
     | B(4,0) | B(4,1) | B(4,2) | B(4,3) | B(4,4) |
     +--------+--------+--------+--------+--------+
     | U(3,0) | U(3,1) | U(3,2) | U(3,3) | U(3,4) |
     | B(3,0) |        |        |        | B(3,4) |
     +--------+--------+--------+--------+--------+
     | U(2,0) | U(2,1) | U(2,2) | U(1,3) | U(2,4) |
     | B(2,0) |        |        |        | B(2,4) |
     +--------+--------+--------+--------+--------+
     | U(1,0) | U(1,1) | U(1,2) | U(1,3) | U(1,4) |
     | B(1,0) |        |        |        | B(1,4) |
     +--------+--------+--------+--------+--------+
     | U(0,0) | U(0,1) | U(0,2) | U(0,3) | U(0,4) |
     | B(0,0) | B(0,1) | B(0,2) | B(0,3) | B(0,4) |
     +--------+--------+--------+--------+--------+

The 9 equations for the 9 unknown interior cells are

 U(1,0) + U(1,2) + U(0,1) + U(2,1) -4U(1,1) = h*h*C(1,1)
 U(1,1) + U(1,3) + U(0,2) + U(2,2) -4U(1,2) = h*h*C(1,2)
 U(1,2) + U(1,4) + U(0,3) + U(2,3) -4U(1,3) = h*h*C(1,3)
 U(2,0) + U(2,2) + U(1,1) + U(3,1) -4U(2,1) = h*h*C(2,1)
 U(2,1) + U(2,3) + U(1,2) + U(3,2) -4U(2,2) = h*h*C(2,2)
 U(2,2) + U(2,4) + U(1,3) + U(3,3) -4U(2,3) = h*h*C(2,3)
 U(3,0) + U(3,2) + U(2,1) + U(4,1) -4U(3,1) = h*h*C(3,1)
 U(3,1) + U(3,3) + U(2,2) + U(4,2) -4U(3,2) = h*h*C(3,2)
 U(3,2) + U(3,4) + U(2,3) + U(4,3) -4U(3,3) = h*h*C(3,3)
 
assigning the 9 unknowns new subscripts (i-1)+3*(j-1) the simultaneous equations,
in matrix form, become:

    (0) (1) (2) (3) (4) (5) (6) (7) (8)   (9)final   (9) initial
     11  12  13  21  22  23  31  32  33    unknown    constant 
   +------------------------------------+ +------+ +------------------------+
11 | -4   1   0   1   0   0   0   0   0 | |U(1,1)| |h*h*C(1,1)-B(1,0)-B(0,1)|
12 |  1  -4   1   0   1   0   0   0   0 | |U(2,1)| |h*h*C(1,2)-B(0,2)       |
13 |  0   1  -4   0   0   0   0   0   0 | |U(3,1)| |h*h*C(1,3)-B(1,4)-B(0,3)|
21 |  1   0   0  -4   1   0   1   0   0 | |U(2,1)| |h*h*C(2,1)-B(2,0)       |
22 |  0   1   0   1  -4   1   0   1   0 |*|U(2,2)|=|h*h*C(2,2)              |
23 |  0   0   1   0   1  -4   0   0   1 | |U(2,3)| |h*h*C(2,3)-B(2,4)       |
31 |  0   0   0   1   0   0  -4   1   0 | |U(3,1)| |h*h*C(3,1)-B(3,0)-B(4,1)|
32 |  0   0   0   0   1   0   1  -4   1 | |U(3,2)| |h*h*C(3,2)-B(4,2)       |
33 |  0   0   0   0   0   1   0   1  -4 | |U(3,3)| |h*h*C(3,3)-B(4,3)-B(3,4)|
   +------------------------------------+ +------+ +------------------------+
   
Notice that the diagonal elements are all  -4.
In a larger matrix it would be more obvious that there are four bands
of  1's  with the remainder of the matrix 0's. Typically a specialized
sparse matrix or band matrix equation solver would be used to solve large
systems of equations.

The "C" language code to generate the ut matrix (nx=5, ny=5) is

static void init_matrix()
{
  int i, j, ii, jj, k;
  
  k = (nx-2)*(ny-2); /* constant row subscript */
  
  /* zero internal cells */
  for(i=1; i<nx-1; i++)
    for(j=1; j<ny-1; j++)
      for(ii=1; ii<nx-1; ii++)
        for(jj=1; jj<ny-1; jj++)
          ut[(i-1)+(ny-2)*(j-1)][(ii-1)+(ny-2)*(jj-1)] = 0.0;
          
  for(i=1; i<nx-1; i++) /* inside boundary */
  {
    for(j=1; j<ny-1; j++)
    {
      ii = (i-1)+(ny-2)*(j-1); /* row of matrix */
      ut[ii][ii] = -4.0;
      ut[ii][k] = h*h*C(xmin+i*h, ymin+j*h);
      if((j+1)<(ny-1))
         ut[ii][(i-1)+(ny-2)*((j-1)+1)] = 1.0;
      else
         ut[ii][k] = ut[ii][k] - B(xmin+i*h, ymin+(j+1)*h);
      if((j-1)!=0)
         ut[ii][(i-1)+(ny-2)*((j-1)-1)] = 1.0;
      else
         ut[ii][k] = ut[ii][k] - B(xmin+i*h, ymin+(j-1)*h);
      if((i+1)<(nx-1))
         ut[ii][((i-1)+1)+(ny-2)*(j-1)] = 1.0;
      else
         ut[ii][k] = ut[ii][k] - B(xmin+(i+1)*h, ymin+j*h);
      if((i-1)!=0)
         ut[ii][((i-1)-1)+(ny-2)*(j-1)] = 1.0;
      else
         ut[ii][k] = ut[ii][k] - B(xmin+(i-1)*h, ymin+j*h);
    }
  }
}

Solving the system of linear equations yields the numerical solution
of the PDE at the grid points  i,j  x=xmin+i*h  y=ymin+j*h