function U = hs2d(F,k,Dx,Dy)
%HPS2D  2D Helmholtz equation solver.
%  U = FHS2D(F,K,DX,DY) solves the Helmholtz equation in two-space 
%  dimensions with source term F(X,Y) defined on a rectangular domain 
%  assuming Dirichlet boundary conditions and for constant wavenumber K.  
%  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 PS2D, FPS2D, FHS2D.

%  Francisco J. Beron-Vera, 2006/11/22

% 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 + k^2*U = F by Gaussian elimination.
k2 = spdiags(repmat(k^2,Nx*Ny,1),0,Nx*Ny,Nx*Ny);
U(2:end-1,2:end-1) = reshape((nab2+k2)\F,Ny,Nx)