function [x, y, u, residual, iter] = poisson(n,f,g,tol,maxiter)
% poisson Numerically approximates the solution of the Poisson
% equation on the square 0 <= x,y <= 1 using the preconditioned
% conjugate method.
%
%   [x y u residual iter] = poisson(n,f,g,tol,numiter) computes
%   the solution. n is the number of subintervals, f is the right-hand
%   side, g is the boundary condition on the square, tol is the required
%   error tolerance, and maxiter is the maximum number of iterations
%   to execute. x and y are the grid of points that form the
%   approximate solution in the interior of the square, and u is
%   the matrix containing the numerical solution. residual reflects
%   the accuracy of the computation obtained after application of precg.
%   if the iteration succeeds in ataining the error tolerance, iter
%   is the number of iterations needed, or -1 if the iteration failed
%   to attain the tolerance. if the iteration succeeded, poisson draws
%   a surface plot of the solution; otherwise it outputs an error message.

h = 1/n;
m = (n-1)^2;
u = zeros(n-1,n-1);
x = (h:h:1-h)';
y = x;
P = gallery('poisson',n-1);
b = zeros(m,1);
k = 1;
% initialize vector b with values of f
for i = 1:n-1
    for j = 1:n-1
        b(k) = h^2*feval(f,x(i),y(j));
        k = k + 1;
    end
end
% now add the boundary condition u = g(x,y) on boundary of the square
b(1) = b(1) + feval(g,0,h) + feval(g,h,0);
% fprintf('Assign b(1) = b(1) + g(0,h) + g(h,0)\n');
b(n-1) = b(n-1) + feval(g,1-h,0) + feval(g,1,h);
% fprintf('Assign b(%d) = b(%d) + g(1-h,0) + g(1,h)\n',n-1,n-1);
b(m-(n-2)) = b(m-(n-2)) + feval(g,0,1-h) + feval(g,h,1);
% fprintf('Assign b(%d) = b(%d) + g(0,1-h) + g(h,1)\n',m-(n-2),m-(n-2));
b(m) = b(m) + feval(g,1-h,1) + feval(g,1,1-h);
% fprintf('Assign b(%d) = b(%d) + g(1-h,1) + g(1,1-h)\n',m,m);

% bottom interior
j = 1;
for i = 2:n-2
    b(i) = b(i) + feval(g,j*h,0);
    j = j + 1;
%     fprintf('Assign b(%d) = b(%d) + g(%f,0)\n',i,i,i*h);
end

% top interior
j = 1;
for i = m-(n-3):m-1
    b(i) = b(i) + feval(g,j*h,1);
    j = j + 1;
%     fprintf('Assign b(%d) = b(%d) + g(%f,1)\n',i,i,i*h);
end

% left interior
j = 1;
for i = n:n-1:m-(n-2)-1
    b(i) = b(i) + feval(g,0,j*h);
    j = j + 1;
%     fprintf('Assign b(%d) = b(%d) + g(0,%f)\n',i,i,i*h);
end

% right interior
j = 1;
for i = 2*(n-1):n-1:m-n
    b(i) = b(i) + feval(g,1,j*h);
    j = j + 1;
%     fprintf('Assign b(%d) = b(%d) + g(1,%f)\n',i,i,i*h);
end

[z, residual, iter] = precg(P,b,zeros(length(b),1),tol,maxiter);
k = 1;
for i = 1:n-1
    for j = 1:n-1
        u(i,j) = z(k);
        k = k + 1;
    end
end

if iter ~= -1
    % Plot the result
    [X, Y] = meshgrid(x,y);
    surf(X,Y,u);
else
    disp('Not sufficient accuracy for a plot');
end

