Dirichlet problem for 2D Poisson equation - finite difference method

A simple example
- regular mesh on a square, homogenneous boundary condition, right hand side of the equation identically equal to 1 (or 1 on one half of the square and -1 on the other half, respectively):
n=20;             % n - number of steps at directions x, y
G=numgrid('N',n); % G - the mesh
D=delsq(G);       % D - the system matrix
N=sum(G(:)>0);    % N - the number of unknowns
F=ones(N,1);      % F - the right hand side
% F(1:round(N/2))=-1; % an alternative - right hand side changed
u=D\F;            % u - the solution
  % graph of the solution:
V=G;
V(G>0)=full(u(G(G>0)));
mesh(V)
you should obtain this, or that figure, respectively.

Look at the system matrix (for small n):
n=5;              % n - number of steps at directions x, y
G=numgrid('N',n); % G - the mesh
D=delsq(G);       % D - the system matrix
full(D)           % shows the matrix D on the screen