function u = fft_poisson(b,h)
% This function solves the 2-d Poisson problem del^2 u = f(x,y) using the fast
% Fourier transformation. The Poisson problem has homogeneous Dirichlet boundary
% conditions, and is defined on a square region with equal sized mesh widths.
%
% parameters:
% b = matrix of f values evaluated at interior meshpoints
% h = mesh width
%
% returns:
% u = solution to PDE at interior meshpoints
% get dimensions of b
[m, n] = size(b);
b_bar = zeros(n,n);
u_bar = b_bar;
u = u_bar;
% make sure we have a square grid
if (m ~= n)
   error('matrix b is not square');
end
% b_bar = 2/(n+1) * v * b * v
% first do a DST on columns of b, which is the same as multiplying v*b
for j=1:n
  b_bar(:,j) = fast_dst(b(:,j));
end
% then do DST on rows of vb, which is analogous to multiplying v*b*v
for i=1:n
  % have to take transpose of row, since fast_dst needs a column vector
  b_bar(i,:) = fast_dst(b_bar(i,:)')';
end
% now scale by 2/(n+1)
b_bar = b_bar * (2/n+1);
% next we can solve for u_bar
u_bar = zeros(n,n);
lambda = [1:n];
lambda = -4 * (sin((lambda*pi) / (2*n + 2))).^2;
for i=1:n
 for j=1:n
  u_bar(i,j) = (h^2 * b_bar(i,j)) / (lambda(i) + lambda(j));
 end
end
% u = 2/(n+1) * v * u_bar * v
% do a DST on columns of u_bar, which is analogous to multiplying v * u_bar
for j=1:n
 u(:,j) = fast_dst(u_bar(:,j));
end
% then do a DST on rows
for i=1:n
  u(i,:) = fast_dst(u(i,:)')';
end
% last, multiply by 2/(n+1)
  u = u * 2/(n+1);
return