function U = fps2d(F,Dx,Dy)
%FPS2D  Fast 2D Poisson equation solver.
%  U = FPS2D(F,DX,DY) solves the Poisson equation in two-space dimensions 
%  with source term F(X,Y) defined on a rectangular domain and Dirichlet 
%  boundary conditions. The boundary terms are specified in the perimeter
%  of the NY x NX matrix F generated using meshgrid. Thus EAST, WEST,
%  SOUTH, and NORTH boundary conditions are to be specified, respectively, 
%  in F(:,1) F(:,END), F(1,:), and F(END,:). Meshgrid widths in X and Y 
%  directions are specified in DX and DY, respectively. Inversion of the 
%  discrete five-point Laplacian is perfomed using the fast sine transform.
%  For best performance speed NX - 1 and NY - 1 should be chosen as powers 
%  of 2.
%
%  See also PS2D, HS2D, FHS2D.

%  Francisco J. Beron-Vera, 2006/11/15
%  $Revision: 1.1 $  $Date: 2006/11/22 21:11:20 $

% Get dimensions.
Nx = size(F,2)-2;
Ny = size(F,1)-2;

% Fold boundary into source term.
U = F;
U(2:end-1,2:end-1) = 0;
U(2:end-1,2:end-1) = 0;
F = F - 4*del2(U,Dx,Dy);
F([1 end],:) = [];
F(:,[1 end]) = [];

% Eigenvalues of discrete Laplacian (nab2).
[nx ny] = meshgrid(1:Nx,1:Ny);
L = 2/Dx^2*(cos(nx*pi/(Nx+1))-1) + 2/Dy^2*(cos(ny*pi/(Ny+1))-1);

% Solve nab2*U = F using fast sine transform.
U(2:end-1,2:end-1) = idst(idst(dst(dst(F)')'./L)')';