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