%
%  code01.m -- linear dispersion
%
%		u  + d u    = 0
%		 t      xxx
%
%               IVs:  u(x,0) = f(x)
%               BCs:  u -> 0 as x -> +/- infinity
%
%  fourier inversion (period L)
%
%  u(x,t) = ifft[ fft[f(x)] e^(i k^3 t) ]
%
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

%  two cases: choice = 1 or 2

choice = 1;

switch choice
  case {1}
    d =  1      ;
    L = 16*pi   ;  N  = 512   ;
    Tend = 0.25 ;  dt = Tend/4;
    blurb = '\bf dispersion:  u_t + u_{xxx} = 0';
  case {2}
    d = -1      ;
    L = 16*pi   ;  N  = 512   ;
    Tend = 0.25 ;  dt = Tend/4;
    blurb = '\bf dispersion:  u_t - u_{xxx} = 0';
end

%  set grids: wavenumbers & coefficients
k  = (2*pi/L)*(-N/2:N/2-1)';

dx = L/N;  xg = (0:dx:L-dx)';  x = [xg ; L];

%  IV initialization & FFT
switch choice
  case {1}
    IVs  = exp(-2*(xg-(L-pi)).^2);  pmin = -0.5;  pmax = 1.2;
  case {2}
    IVs  = exp(-2*(xg-   pi ).^2);  pmin = -0.5;  pmax = 1.2;
end

t = 0;  ug = IVs;  u = [ug ; ug(1)];
clf;  hold on;  axis([0 L pmin pmax]);
plot(x,u,'r');  plot(x,u,'k.','markersize',7);

vg   = fftshift(fft(ug))/N;

%  fourier inverse solution

for time = 0:dt:Tend
  fact = exp(i*d*time*(k.^3));
  ug = real(ifft(fftshift(vg.*fact))*N);
  plot(x,[ug ; ug(1)],'b');
  pause(0.1)
%  pause
end

plot(x,[ug ; ug(1)],'r');  
plot(x,[ug ; ug(1)],'k.','markersize',7);
title(blurb);  xlabel('\bf x-axis');  ylabel('\bf u-axis')

disp('hit return for blow-up')
pause

switch choice
  case {1}
    axis([0.75*L      L pmin pmax])
    text(L-pi,1.1,'\bf gaussian IV')
  case {2}
    axis([   0*L 0.25*L pmin pmax])
    text(  pi,1.1,'\bf gaussian IV')
end