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

% defining mesh size and grid coordinates

  h = .2;                         % mesh size
  x = 0:h:3;                      % h separated x coordinates from 0 to 3
  y = -1:h:1;                     % corresponding y coordantes
  [X,Y] = meshgrid(x,y);          % extend to matrices of x and y coordinates


% defining boundary values

  g = cos(Y.*(2-X));              % desired boundary values (note the .*)


% solving using "Svenssons (5-point) formula"

  u = g;                          % initializing u values
  I = 2:size(u,1)-1;              % first indices for interior mesh points
  J = 2:size(u,2)-1;              % corresponding second indices
  u(I,J) = 0;                     % initializing to zero in the interior 
  res = 1;                        % initializing "residual"
  set(gcf,'Doublebuffer','on')    % for smooth plot updating

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

    % update (all) interior values according to Svenssons rule

    u(I,J) = (u(I+1,J) + u(I-1,J) + u(I,J+1) + u(I,J-1)) / 4;

    % compute the residual error (seek to get res = 0)

    res = u(I+1,J) + u(I-1,J) + u(I,J+1) + u(I,J-1) - 4 * u(I,J);

    surf(X,Y,u)                   % plot current u
    axis equal
    xlabel('x')
    ylabel('y')
    drawnow
    shg                           % bring graph window up front

  end