  function U = BFastPoisson(a,b,n1,c,d,n2,f)
% Solves the Poisson problem
%
%                u_xx(x,y) + u_{yy}(x,y) = 0
%
% over the rectangle R = {(x,y) | a<=x<=b, c<=y<=d } with
% u(x,y) = f(x,y) on the boundary of R.
%
% U is (n1+1)x(n2+1) where U(i,j) \approx u(x(j),y(i)) where 
%
%           x(j) = a + (j-1)hx   , j=1:n2+1 where   
hx = (b-a)/(n2+1);
%           y(i) = d + (i-1)hy   , i=1:n1+1 where
hy = (c-d)/(n1+1);
% Note that x = linspace(a,b,n2+1) and y = linspace(d,c,n1+1)


U = zeros(n1+1,n2+1);
% The boundary values...
for j=1:n2+1
    U(1,j)    = f(a+(j-1)*hx,d);
    U(n1+1,j) = f(a+(j-1)*hx,c);
end
for i=1:n1+1
    U(i,1)    = f(a,d-(i-1)*hy);
    U(i,n2+1) = f(b,d-(i-1)*hy);
end
fMatrix = zeros(n1-1,n2-1);
% At each mesh point we have (in compass notation):
%     4u(P) - u(North) - u(East) - u(South) - u(West) = 0
% Along the "edges" some of P's neighbors are on the boundary and these
% known values move over to the right hand side.

fMatrix(1,:)    = U(1,2:n2);
fMatrix(n1-1,:) = U(n1+1,2:n2);
fMatrix(:,1)    = fMatrix(:,1)    + U(2:n1,1);
fMatrix(:,n2-1) = fMatrix(:,n2-1) + U(2:n1,n2+1);

U(2:n1,2:n2) = FastPoissonSolve(fMatrix,hx,hy);