function [soln, xs, ys] = poissonSolver(a, b, h, ff, lBC, rBC, bBC, tBC)
%
% function [soln, xs, ys] = poissonSolver(a, b, h, ff, lBC, rBC, bBC, tBC)
%
% written by Thomas L. Scofield on Nov. 17, 2009
%
% This function uses standard finite differences and a uniform
% mesh to solve Poisson's equation
%     - Laplacian u = f
% subject to Dirichlet boundary conditions in a rectangular region,
% 0 < x < a, 0 < y < b, of the plane:
%     u(0,y) = lBC(x)
%     u(a,y) = rBC(x)
%     u(x,0) = bBC(x)
%     u(x,b) = tBC(x)
%
% INPUTS:
%    a, b  specify top max x and y-values for the rectangle
%          in which the problem is posed.
%    h     width between mesh points in both x and y-directions
%    ff    forcing function
%    lBC   a function which provides values along the left boundary x = 0
%    rBC   a function which provides values along the right boundary x = a
%    bBC   a function which provides values along the bottom boundary y = 0
%    tBC   a function which provides values along the top boundary y = b
%
% OUTPUT:
%    soln  a matrix of solution values, corresponding to inputs xs and ys
%    xs    a matrix of x-values
%    ys    a matrix of y-values
%          those resulting from the command [xs, ys] = meshgrid(0:h:a, 0:h:b);

  enn = a / h - 1;
  emm = b / h - 1;
  [xs, ys] = meshgrid(0:h:a, 0:h:b);

  soln = zeros(emm+2, enn+2);
  % incorporate boundary data into the solution
  soln(1:emm+2, 1) = lBC(ys(1:emm+2,1));
  soln(1:emm+2, enn+2) = rBC(ys(1:emm+2,1));
  soln(1, 1:enn+2) = bBC(xs(1,1:enn+2))';
  soln(emm+2, 1:enn+2) = tBC(xs(1,1:enn+2))';

  % set up right-hand side vector
  dat = h^2 * ff(xs(2:emm+1,2:enn+1), ys(2:emm+1,2:enn+1));
  dat = dat';
  dat = dat(:);
  dat(1:enn) = dat(1:enn) + soln(1, 2:enn+1)';
  dat((emm - 1)*enn + (1:enn)) = dat((emm - 1)*enn + (1:enn)) ...
                                  + soln(emm+2, 2:enn+1)';
  dat(1:enn:emm*enn) = dat(1:enn:emm*enn) + soln(2:emm+1, 1);
  dat(enn:enn:emm*enn) = dat(enn:enn:emm*enn) + soln(2:emm+1, enn+2);

  A = discreteLaplacian(enn, emm);
% A = badDiscreteLaplacian(enn, emm);

% B = full(A);
  savetime = cputime();
  uu = A \ dat;
  cputime() - savetime
% savetime = cputime();
% uu = B \ dat;
% cputime() - savetime
  for j = 2:emm+1
    soln(j, 2:enn+1) = uu((j - 2)*enn + (1:enn));
  end
  mesh(xs, ys, soln)
end