function X = FastPoisson( n, f )
% Fast Poisson solver using 2nd order central differences.
% 
%  Solves    -laplacian(u) = f   with zero Dirichlet on 
%               on [0 1]x[0 1]   with an nxn grid. 
%
% Alex Townsend, September 2014.

% Evaluate rhs on gird:
x = linspace(0, 1, n+2); x([1,end])=[];
[xx, yy] = meshgrid(x);
F = f(xx, yy) ;

% Set up and solve matrix equation KX + XK = F:
h = 1/n; s = sqrt(2/(n+1));
d = (2-2*cos((1:n)'*pi/n))/h^2;     % Eigenvalues of K.
F = s*dst1(s*dst1(F).').';          % S \ F * S
Y = F./((d-1000)*ones(1,n) + ones(n,1)*d');% Divide by eigvals
X = s*dst1(s*dst1(Y).').';          % X = S \ Y * S

% Add back in zero dirichlet conditions
Z = zeros(1,n);
X = [0   Z   0 ;
     Z'  X   Z';
     0   Z   0 ];
X = real( X );

% Plot
x = linspace(0, 1, n+2);
[xx, yy] = meshgrid(x);
surf(xx, yy, X,'edgealpha',0,'facecolor', 'interp'),
view(0,90), axis square, set(gca, 'fontsize', 16)
end