m=15;                            % number of gridpoints
A=Laplacian(m,2);                % five point Laplacian
f=zeros(m*m,1); f(m:m:end)=1;    % put bc into the rhs
u=A\f;                           % solve by sparse Gaussian elimination
x=0:1/(m+1):1; y=x;              % mesh point vectors
ugs=zeros(size(f));              % initial guess
U=zeros(m+2); U(end,1:m+2)=1;    % for plotting
LD=tril(A);                      % L+D
for k=0:10                       % do 10 Gauss-Seidel steps
  err(k+1)=max(max(abs(u-ugs))); % compute error				
  U(2:m+1,2:m+1)=reshape(ugs,m,m);  
  mesh(x,y,U);                   % plot current approximation
  xlabel('x');ylabel('y');
  ugs=ugs+LD\(f-A*ugs);          % perform one Gauss-Seidel step
  pause
end
