function [ xi, yi, u ] = Poisson2D( n, f )
%POISSON2D FFT method for solving 2D Poisson equation on a square
%   Using standard O(h^2) finite differences we get a linear
%   (block circulant with circulant blocks) system of equations
%   to solve for u: -(u_xx + u_yy) = f gives A u = f.
%   Using the 2d FFT diagonalizes the n^2 by n^2 system.

% Use left/upper part of grid
x = (0:n)/n;
h = 1/n;
[xi,yi] = meshgrid( x );
fv = feval(f, xi, yi );
% zero out the outer edges
fv(1,:) = zeros(1,n+1);
fv(n+1,:) = fv(1,:);
fv(:,1) = fv(1,:)';
fv(:,n+1) = fv(:,1);
% Odd extension
fv = [fv, -fv(:,end-1:-1:2); -fv(end-1:-1:2,:), fv(end-1:-1:2,end-1:-1:2) ];
% Henrici's Fn is Matlab's ifft
fhat = ifft2( fv );
% Explicit computation of dhat by formula
[k, m] = meshgrid( 2*(n):-1:1 );
dhat = (1-cos(k*pi/n)/2 - cos(m*pi/n)/2)*h^2;
dhat(1,1) = 1; % Don't divide by zero.  
fz = fhat./dhat; % Nans will show up otherwise, and they hurt.
u = h^2*fft2( fz )/4/n^2;
u = real(u(1:n+1,1:n+1)); % Crude code, real u only
% Explicitly impose BC, avoid roundoff error
u(1,:) = zeros(1,n+1);
u(n+1,:) = u(1,:);
u(:,1) = u(1,:)';
u(:,n+1) = u(1,:)';

end

