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