function help_Poisson

N=20;								% Number of points 
L=1;
h=L/N;                              % Grid spacing

X=linspace(0,L,N+1);                   % x-coordinates length N+1 at grid nodes
Y=linspace(0,L,N+1);                   % y-coordinates length N+1 at grid nodes

% Assemble x and y - coordinate at cell centers into 1-D vector length N*N

[x,y]=meshgrid((X(1:end-1)+X(2:end))/2,(Y(1:end-1)+Y(2:end))/2);
x=reshape(x,[],1);
y=reshape(y,[],1);

% Define the forcing function

f=     ;                         % Forcing function


% Set up matrix A (including boundary condition) 
% to solve the system of equations directly 

A=

% Solve the equation A u = h^2 f
% Hint: use \ (backslash) operator (see Matlab documentation: help mldivide)

p_exact=

% Set up u for iterative mthod

p=     ;                         % Inital condition


% Relaxation parameter

omega= ;                         % =1 => Gauss-Seidel   =1.5 =>  S.O.R.

% Choose a criterion for convergence and the maximum iteration steps

err=1e-8;
Mmax=5000;
m=0;

% March trough all the inner points
% p(X_i,Y_j)=(1-omega)*p(X_i,Y_j)+omega/4*(p(X_i,Y_(j-1))+p(X_i,Y_(j+1))+p(X_(i-1),Y_j)+p(X_(i+1),Y_j)+h^2*f(X_i,Y_j))
% Remember that p now should be a 1-D vector 

while convergence>err &&  m<Mmax 
    m=m+1;              

	% i=1, j=1
	k=1;
	p(k)=

	% along i=1 
	for j=2:(N-1)                      
		k=j;
		p(k)=
	end

	% i=1, j=N
	k=N;
	p(k)=

    for i=2:(N-1)                      
		% j=1
		k=1+N*(i-1);
		p(k)=

        for j=2:(N-1)
            k=j+N*(i-1);
            p(k)=
        end

		% j=N
		k=N+N*(i-1);
		p(k)=

    end
    
	% i=N, j=1
	k=1+N*(N-1);
	p(k)=

	% along i=N
	for j=2:(N-1)                       
		k=j+N*(N-1);
		p(k)=
	end

	% i=N, j=N
	k=N*N;
	p(k)=

    % Calculate Residuals
    
    R(m)=
    convergence=R(m);
    
    % Calculate difference to exact solution
    
    epsilon(m)=
    
end


surf(X,Y,reshape(p,N,N))
xlabel('X'),ylabel('Y')
figure(2)
semilogy(R);

