1 % laplace.m
 2 % Author: Lorenzo Bolla <lbolla (at) gmail dot com>
 3 % Last modified: Apr. 11, 2007
 4  
 5 xmin = 0;
 6 xmax = 1;
 7 ymin = 0;
 8 ymax = 1;
 9  
 10 nx = 200;
 11 ny = nx;
 12  
 13 dx = (xmax - xmin) / (nx - 1);
 14 dy = (ymax - ymin) / (ny - 1);
 15  
 16 x = xmin:dx:xmax;
 17 y = ymin:dy:ymax;
 18  
 19 n_iter = 100;
 20 eps = 1e-16;
 21  
 22 u = zeros(nx,ny);
 23 u(1,:) = xmin^2 - y.^2;
 24 u(end,:) = xmax^2 - y.^2;
 25 u(:,1) = x.'.^2 - ymin^2;
 26 u(:,end) = x.'.^2 - ymax^2;
 27 err = Inf;
 28 count = 0;
 29 while err > eps && count < n_iter
 30  count = count + 1;
 31  dx2 = dx^2;
 32  dy2 = dy^2;
 33  dnr_inv = 0.5 / (dx2 + dy2); 
 34  u_old = u;
 35  u(2:end-1, 2:end-1) = ((u(1:end-2, 2:end-1) + u(3:end, 2:end-1))*dy2 + (u(2:end-1,1:end-2) + u(2:end-1, 3:end))*dx2)*dnr_inv;
 36  v = (u - u_old);
 37  err = sqrt(v(:).'*v(:));
 38 end