% Script for solving u'' = exp(4x) , u(0) = u(1) = 0 using spectral method
% based on N-point DFT using odd extension to 0 < x < 2. Solution is plotted
% for N = 4 (h = 0.5) and max-norm error is plotted vs. h = 2/N for
% N = 4, 8, 16, 32, 64

  L = 4;
  xl = -1;
  xr = 3;
  pmax = 5;
  for p = 1:pmax
    N = 2^(p+1);
    h(p) = L/N;
    k = [0:(N/2-1) (-N/2):(-1)]*2*pi/L;
    k(1) = 1e-6;          % Redefine k(1) nonzero to prevent 0/0 division later

    x = xl + (0:(N/2))*h(p);    % These are only the grid points covering [0 1]
     
%    Define right hand side over extended domain 
    
    f = zeros(1,N);
    f(1) = 0;
    f(N/2+1) = 0;
    jrange = 2:(N/2);
    f(jrange) = exp(4*x(jrange));
    f(N + 2 - jrange) = -f(jrange);       %
       
%    DFT, solve, and inverse DFT

    fhat = fft(f);
    u = real(ifft(-fhat./k.^2));

% For N = 16 (medium resolution) save spectral solution for plotting

    if p == 3
      u1 = u;
    end

% Restrict u to gridpoints x in [-1 1]

    u = u(1:(N/2+1));            

%    Calculate error 

    uex = (exp(4*x)-x*sinh(4) -cosh(4))/16;
    err(p) = norm(u-uex,inf);
  end

%  Plot N = 16 solution using interpft to interpolate to our high-resolution
%  (N = 64) grid.

  subplot(2,2,1)
  loglog(h,err,'x',h,err(5)*(h/h(5)).^2)
  xlabel('h')
  ylabel('Maximum error')
  title('Spectral method on u" = exp(4x)')
  legend('Error','h^2 fit',2)

  subplot(2,2,2)
  u1hi = interpft(u1,N);
  plot(x,uex,'-',x,u1hi(1:(N/2+1)),'-.');
  xlabel('x')
  ylabel('u')
  title('N = 16')
  legend('Exact','Spectral',0)