%
% solves poisson equation in [0,1]x[0,1] using (n+2)x(n+2) grid.
% using sor method.  1<=ome<2
%
function poisson_dir(n,ome)
h=1/(n+1);
h2=h*h;
xx=0:h:1;
yy=0:h:1;
[X,Y]=meshgrid(xx,yy);
f=X.*Y;
%f=ones(n+2,n+2);
u=zeros(n+2,n+2);
itmax=2000;
for iter=1:itmax
    uold = u;
    for j=2:n+1
        for i=2:n+1
            d = 4.0; 
            uu = (h2*f(i,j)+u(i+1,j)+u(i-1,j)+u(i,j+1)+u(i,j-1))/d;
            u(i,j) = (1-ome)*u(i,j) + ome*uu;
        end
    end
    if norm(u-uold,'inf') < 1e-6*norm(u,'inf')
        break;
    end
end
iter
max(max(u))
contourf(X,Y,u,15); axis equal; colorbar