function solution = SimplePoisson(xVals,yVals,forcing)

%
% Function that finds the solution to the two-dimensional
% Poisson equation on a rectangular grid when the
% interior "forcing" and solution on the boundaries is
% given.
%
%     solution = SimplePoisson(xVals,yVals,forcing)
%
% Input:
%
%   xVals: A one dimensional array of size n containing
%          the values of the first independent variable
%          on the grid points.  The values must be
%          monotonously increasing.
%
%   yVals: A one dimensional array of size m containing
%          the values of the second independent variable
%          on the grid points.  The values must be
%          monotonously increasing.
%
% forcing: An m x n array that contains info about the
%          desired solution.
%          * The value of forcing(i,j) at any interior
%          point of the grid should be the "forcing" at
%          that point (i.e. whatever it is at the point
%          that wants to make the solution nontrivial).
%          * The value of forcing(i,j) at any boundary
%          point, (i=1, i=m, j=1, or j=n), should be the
%          desired value of the solution at that point.
%
% Output:
%
% solution: An m x n array containing the values of the
%           approximate solution at the m x n points.
%
% Hint: arrays forcing and solution must be compatible
% with the grid points that you would get from
%
%      [xGrid yGrid] = meshgrid(xVals,yVals)
%
% In particular an array element with indices (i,j)
% corresponds physically to position yVals(i),xVals(j).
% In still other words, when the array is displayed, x
% increases from left to right and y from top to bottom.

% check array xVals and find n
if ~isvector(xVals)
    fprintf('*** SimplePoisson: ')
    error('Array xVals is not one-dimensional!')
end
n = max(size(xVals))
% check array yVals and find m
if ~isvector(yVals)
    fprintf('*** SimplePoisson: ')
    error('Array yVals is not one-dimensional!')
end
m = max(size(yVals))
% check array forcing
if ~isequal(size(forcing),[m n])
    fprintf('*** SimplePoisson: ')
    error('Array forcing has the wrong size!')
end

% check that we have interior points
if m<3 | n<3
    fprintf('*** SimplePoisson: ')
    error('There are no interior points!')
end

% compute the number of unknown solution values
N=m*n
% initialize storage for the diagonals of the matrix
Diags=zeros(N,5);
Diags(:,3)=ones(N,1);
% zero the vector of right hand sides
b=zeros(N,1);

% form the interior equations for all y positions
for i = 2:m-1
    % compute the mesh increments in y
    dy1=yVals(i)-yVals(i-1);
    dy2=yVals(i+1)-yVals(i);
    if dy1<=0 || dy2<=0
        fprintf('*** SimplePoisson: ')
        error('The yVals values are not increasing!')
    end
    % find the average increment
    dyAve=.5*(dy1+dy2);
    % form the equations for all x positions on the line
    for j = 2:n-1
        % compute the equation number
        I = i + (j-1)*m;
        % compute the mesh increments in x
        dx1=xVals(j)-xVals(j-1);
        dx2=xVals(j+1)-xVals(j);
        if dx1<=0 || dx2<=0
            fprintf('*** SimplePoisson: ')
            error('The xVals values are not increasing!')
        end
        % find the average increment
        dxAve=.5*(dx1+dx2);
        % diagonal coefficient of the equation
        Diags(I,3)=-1/(dx2*dxAve)-1/(dx1*dxAve)...
                   -1/(dy2*dyAve)-1/(dy1*dyAve);
        % off-diagonal coefficients (stupid Matlab order)
        Diags(I+m,5)=1/(dx2*dxAve);
        Diags(I-m,1)=1/(dx1*dxAve);
        Diags(I+1,4)=1/(dy2*dyAve);
        Diags(I-1,2)=1/(dy1*dyAve);
        % right hand side of the equation
        b(I)=forcing(i,j);
    end
end

% combine the diagonals into a sparse NxN matrix
A=spdiags(Diags,[-m -1 0 1 m],N,N);

% add the x-boundary conditions at i = 1 and m
for j = 1:n
    % x-boundary condition at i = 1
    I = 1 + (j-1)*m;
    %A(I,I)=1;
    b(I)=forcing(1,j);
    % x-boundary condition at i = m
    I = m + (j-1)*m;
    %A(I,I)=1;
    b(I)=forcing(m,j);
end

% add the y-boundary conditions at j = 1 and n
for i = 2:m-1
    % y-boundary condition at j = 1
    I = i;
    %A(I,I)=1;
    b(I)=forcing(i,1);
    % y-boundary condition at j = n
    I = i + (n-1)*m;
    %A(I,I)=1;
    b(I)=forcing(i,n);
end

% check the condition number of the matrix
if N <= 1000
    condA=cond(A);
    solRelerrDueToMatlab=condA*eps
else
    disp('Using an estimated L1 condition number.')
    condA1=condest(A);
    solRelerrDueToMatlab=condA1*eps
end
if solRelerrDueToMatlab >= 0.01
    fprintf('*** SimplePoisson: The equations are ')
    error('singular!  Check xVals and yVals.')
end

% find the solution as a vector
solutionVec=A\b;
% reformat it to the same form as the forcing array
solution=reshape(solutionVec,[m n]);

end
