% Numerically calculate soln. to advection eqn. dq/dt + dq/dx = 0 on
% domain 0<x<1 with periodic BCs using spectral method with RK4
% time-differencing scheme.

  N = 8; % Number of modes
  nu = 0.25; % Courant number 
  L = 1;  % Domain size
  x = L*(0:(N-1))/N;  %  x-gridpoints [1xN]
  dx = L/N;
  M = [0:(N/2-1) (-N/2):(-1)];  
  k = 2*pi*M/L;  % Wavenumbers [1xN].  

  q0 = 1./(4+3*cos(2*pi*x/L));  % Initial condition
  tf = 1; % Final time
  dt = nu*dx; % Timestep
  nt = round(tf/dt); % Number of timesteps to take

  qhat = fft(q0);  % Initial Fourier expansion coeffs [1xN].
  for it = 1:nt

    % March forward dqhat/dt = -Shat using RK4
    % where Shat is DFT of S(q) = dq/dx

    d1 = -dt*1i*k.*qhat; 
    d2 = -dt*1i*k.*(qhat + 0.5*d1);
    d3 = -dt*1i*k.*(qhat + 0.5*d2);
    d4 = -dt*1i*k.*(qhat + d3);
    qhat = qhat + (d1 + 2*d2 + 2*d3 + d4)/6; % New qhat
  end
  q = ifft(qhat); % Numerical q at tf [1xN]

  Nf = 256;   % Number of plotting points
  xf = (0:(Nf-1))*L/Nf;
  qf = interpft(q,Nf);  % Numerical solution on plotting grid
  q0f = 1./(4+3*cos(2*pi*xf/L));  % Initial condition
  plot([xf L],[q0f q0f(1)],'b-',[x L],[q q(1)],'r+',[xf L],[qf qf(1)],'r-') 
  xlabel('x')
  ylabel('q')
  axis([0 1 0 1.5])
  title(['Spectral method, RK4, dt/dx = ',num2str(nu),', t = ' num2str(tf)])
  legend('Exact','Gridpoints','Continuous')