function MySvenssonSolver()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Solves the Laplace equation
%
%       div q = 0,  q = - grad u,
%
%  in part of (!) a given rectangular domain, with boundary conditions
%
%       u = g
%
%  outside
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% defining mesh size and grid coordinates

  h = .1;
  x = 0:h:3;
  y = -1:h:1;
  [X,Y] = meshgrid(x,y);


% defining boundary values

  g = zeros(size(X));                      % desired boundary values                       
  f = ones(size(X));                       % given right hand side


% solving using Svenssons (5-point) formula

  u = g;
  I = 2:size(u,1)-1;
  J = 2:size(u,2)-1;
  u(I,J) = 0;
  res = 1;
  set(gcf,'Doublebuffer','on')

  F = (0 < Y | X > 1);     % setting F = 1 for y > 0 and x > 0, = 0 elsewhere
  F = F(I,J); 

  while max(max(abs(res))) > .002

    u(I,J) = F.*((u(I+1,J) + u(I-1,J) + u(I,J+1) + u(I,J-1)) / 4 ...
           + f(I,J) * h * h / 4);
    res = F .* (f(I,J) * h * h + u(I+1,J) + u(I-1,J) + u(I,J+1) + u(I,J-1) - 4 * u(I,J));

    surf(X,Y,u)
    axis equal
    xlabel('x')
    ylabel('y')
    drawnow
    shg

  end