function MySvenssonSolver()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Solves the Laplace equation
%
%       div q = 0,  q = - grad u,
%
%  in a given rectangular domain, with boundary conditions
%
%       n.q = k (g - u)
%
%  on part of the boundary (and  u = g  on the rest)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 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')

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

    u(I,J) = (u(I+1,J) + u(I-1,J) + u(I,J+1) + u(I,J-1)) / 4 ...
           + f(I,J) * h * h / 4;
  % b.c. n.grad u = k (g - u) with k = 0
    u(I,1) = u(I,2);
  % b.c. n.grad u = k (g - u) with k = 7 and g = y
%    u(I,1) = (u(I,2)/h + 70*Y(I,1)) / (1/h + 70);
    res = 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