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