function [u,x]=FFPoisson1d(f,a,b,n);
% FFPOISSON1D solves the 1d Poisson equation
%   u=FFPoisson1d(f,a,b,n); solves the 1d poisson equation with right
%   hand side function f, boundary values a and b using the fast Fourier
%   transform. Note that f needs to be periodic.

h=(b-a)/n; x=(a:h:b)';                        % finite difference mesh
b=feval(f,x(1:n));
bh=fft(b);                                    % check mean in rhs and
if abs(bh(1))>100*eps,                        % remove if necessary
  warning(['right hand side function mean ' num2str(bh(1)) 'removed']);
  bh(1)=0;
end;
bhint(2:n/2)=-bh(2:n/2)./([1:n/2-1]'.^2);     % ordering of -k^2 in
bhint(n/2+1:n)=-bh(n/2+1:n)./([-n/2:-1]'.^2); % the FFT implementation
bhint(1)=0;                                   % fix average to 0
u=ifft(bhint);
u(n+1)=u(1);                                  % solution is periodic
u=real(u)';                                   % make numerically real
