function U = ps2d(F,Dx,Dy)
%PS2D  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 
%  five-point discrete Laplacian is done by Gaussian elemination as 
%  implemented in Matlab's backslash division.
%
%  See also FPS2D, HS2D, FHS2D.

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

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

% Five-point discrete Laplacian.
DX = repmat(Dx,Nx,1);
DY = repmat(Dy,Ny,1);
d2dx2 = spdiags([1./DX.^2 -2./DX.^2 1./DX.^2],-1:1,Nx,Nx); 
d2dy2 = spdiags([1./DY.^2 -2./DY.^2 1./DY.^2],-1:1,Ny,Ny); 
nab2 = kron(d2dx2,speye(Ny)) + kron(speye(Nx),d2dy2);   

% 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]) = [];
F = F(:);

% Solve nabla2*U = F by Gaussian elimination.
U(2:end-1,2:end-1) = reshape(nab2\F,Ny,Nx)