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