function [x, y, u, residual, iter] = biharmonic(n,f,tol,maxiter)
% BIHARMONIC  Numerically approximates the solution of the biharmonic
% equation on the unit square 0 <= x,y <= 1 using the preconditioned
% conjugate method. The boundary condtions are u(x,y) = 0,
% 0 <= x,y <= 1, normal derivative is zero on the boundary
%
%   [x y u residual iter] = biharmonic(n,f,g,tol,numiter) computes
%   the solution.
% 
%Input:
%  n         - the number of subintervals
%  f         - the right-hand side
%  g         - the boundary condition on the square
%  tol       - the error tolerance
%  maxiter   -  the maximum number of iterations
%
%Output:
%   x, y     - the grid of points that form the approximate solution in the interior of the square
%   u        - the matrix containing the numerical solution
%   residual - the residual
%   iter     - number of iterations required if the computation attains the error tolerance
%                and -1 otherwise
%   plot     - if the iteration succeeds, biharmonic 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);
b = zeros(n-1,n-1);
x = h:h:1-h;
y = x;
for j = 1:n-1
    for i = 1:n-1
        b(i,j) = feval(f,x(i),y(j));
    end
end
b = reshape(b,m,1);
B = biharmonic_op(n-1,1.0,n-1,1.0);
[z, residual, iter] = precg(B,b,zeros(length(b),1),tol,maxiter,...
                            'incomplete Cholesky',1.0e-6);
if iter == -1
    disp('Not sufficient accuracy for a plot');
    return;
end

% form the matrix u from the vector solution z
k = 1;
for i = 2:n
    for j = 2:n
        u(i,j) = z(k);
        k = k + 1;
    end
end

if iter ~= -1
    % Plot the result
   x = 0:h:1;
   y = x;
   [X, Y] = meshgrid(x,y);
   surf(X,Y,u);
   title({'$$\frac{\partial^{4}u}{\partial x^{4}}+2\frac{\partial^{4}u}{\partial x^{2}\partial y^{2}}+\frac{\partial^{4}u}{\partial y^{4}}=f\left(x,\, y\right)\\$$';...
      '$$u=0$, $\frac{\partial u}{\partial n}=0$ on $\partial R=\left\{ 0\leq x,\, y\leq1\right\} $$'},'fontsize', 14, 'FontWeight','bold','interpreter','latex');
   xlabel('x');
   ylabel('y');
   zlabel('z');    
end

