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