function laplace1(N)
%laplace1(N) solves Laplace's equation u_xx + u_yy = 0 on the unit square
%    with (N+1)x(N+1) gridpoints i=1,...,N+1 & j=1,...,N+1
%    BCs are Dirichlet with u=0 on boundary except u(x=1,y)=1
%    Uses Gauss-Seidel iteration

tic
h = 1/N; % dx=dy=h
i = (1:N+1);
x = (i-1)/N;
y = (i-1)/N;
u = zeros(N+1,N+1);
u(i,N+1) = 1; % nonzero BC

sum = 0;
for i = 2:N
    for j = 2:N
        residual = (-4*u(i,j)+u(i+1,j)+u(i-1,j)+u(i,j+1)+u(i,j-1));
        sum = sum+abs(residual);  
    end
end
normresidual0 = sum/(N-1)^2;
normresidual = Inf;
iter = 0;
while normresidual > 10^-5*h^2*normresidual0
    iter = iter + 1;
    sum = 0;
    for i = 2:N
        for j = 2:N
            residual = (-4*u(i,j)+u(i+1,j)+u(i-1,j)+u(i,j+1)+u(i,j-1));
            sum = sum+abs(residual);
            u(i,j) = u(i,j)+residual/4;   
        end
    end
    normresidual = sum/(N-1)^2;
end
iter
toc

figure
surf(x,y,u)
shading flat
xlabel('x','FontSize',16); ylabel('y','FontSize',16); zlabel('u','FontSize',16)

end