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