function U = FastPoisson(a,b,n1,c,d,n2,f)
% Solves the Poisson problem
%
% u_xx(x,y) + u_{yy}(x,y) = -f(x,y)
%
% over the rectangle R = {(x,y) | a<=x<=b, c<=y<=d } with
% u(x,y) = 0 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);
fMatrix = zeros(n1-1,n2-1);
for i=2:n1
for j=2:n2
fMatrix(i-1,j-1) = f(a+(j-1)*hx,d+(i-1)*hy);
end
end
U(2:n1,2:n2) = FastPoissonSolve(fMatrix,hx,hy);