clear all
global N h
global u hhf
global ipts_i ipts_j ipts_n 
N=100;
h=6/(N-3);
etol=0.001^2/N/N/N;

c=1;
for i=1:N
    ipts_n(c)=0;
    for j=1:N
        x=i2x(i);
        y=i2x(j);
        if isint(x,y)
            u(i,j)=0;
            if(ipts_n(c)==0)
                ipts_i(c)=i;
                ipts_j(c)=j;
            end
            ipts_n(c)=ipts_n(c)+1;
        else
            u(i,j)=1;
            if(ipts_n(c)>0)
                c=c+1;
                ipts_n(c)=0;
            end
        end
    end
end
cmax=c-1;
% initialize forcing array
for c=1:cmax
    i=ipts_i(c);
    for j=ipts_j(c):ipts_j(c)+ipts_n(c)-1
        hhf(i,j)=h*h*force(i2x(i),i2x(j));
    end
end

% put the boundary data in
for c=1:cmax
    i=ipts_i(c);
    for j=ipts_j(c):ipts_j(c)+ipts_n(c)-1
        if u(i+1,j)==1
            u(i+1,j)=gbnd(i2x(i+1),i2x(j));
        end
        if u(i-1,j)==1
            u(i-1,j)=gbnd(i2x(i-1),i2x(j));
        end
        if u(i,j+1)==1
            u(i,j+1)=gbnd(i2x(i),i2x(j+1));
        end
        if u(i,j-1)==1
            u(i,j-1)=gbnd(i2x(i),i2x(j-1));
        end
    end
end

% gauss seidel method
for n=1:10000
    uerr=0;
    for c=1:cmax
        i=ipts_i(c);
        for j=ipts_j(c):ipts_j(c)+ipts_n(c)-1
            unew=(u(i-1,j)+u(i+1,j)+...
               u(i,j+1)+u(i,j-1)-hhf(i,j))/4;
            uerr=uerr+(unew-u(i,j))^2;
            u(i,j)=unew;
        end
    end
    if uerr<=etol
        break
    end
end
hold off
pcolor(u);
shading interp 
hold on
contour(u,'k');