% fftpois1d  -u''=f  u(0)=u(1)=0 
clear all
nn=2.^[13:18];
for j=1:size(nn,2)
    n=nn(j); h=1/n; x=[h:h:1-h]'; uex=x.*sin(pi*x);
    tic 
    lam=4*n*n*sin(x*pi/2).^2; f=pi^2*uex-2*pi*cos(pi*x);
    u=imag(fft([0;f;zeros(n,1)]));  %FFT
    u(1)=0; u(n+1:end)=zeros(n,1); u(2:n)=u(2:n)./lam;
    u=imag(fft(u))*2/n; u=u(2:n);  % inverse FFT
    t(j)=toc;
    norm(u-uex,inf),
end
plot(nn.*log(nn),t,'-o'); xlabel('n log n'); ylabel('time');