%
%  lect29.m -- spectral diffusion w/ nonlinearity
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

%  parameters

global n2 N D a

choice = 1

switch choice
  case 1
    D = 0.25;  a = 0.0;
    N =   64;  L =   pi;  dx = 2*L/N;  
    T =    2; dt =  0.4;

  case 2
    D = 0.04;  a =  0.0;
    N =  128;  L =   pi;  dx = 2*L/N;  
    T =    6; dt =   1;

  case 3
    D = 0.04;  a =  0.1;
    N =  128;  L =   pi;  dx = 2*L/N;  
    T =  140; dt =   10;
end

n2 = (( (0:1:N-1)*(pi/2/L) ).^2)';

%  set grids & coefficients

xg = 0:dx:2*L-dx;  x = [xg 2*L];
tg = 0:dt:T;

%  IV initialization (column), FFT & plot

switch choice
  case 1
    IVs  = 1.5*sin(xg)';
  case {2,3}
    IVs  = 1.5*( exp(-4*(xg-2*pi/3).^2) - exp(-9*(xg-7*pi/4).^2 ))';
end

u = IVs;
usave = [u ; u(1)];

u2 = [u ; 0 ; -flipud(u(2:end))];
uf = -imag(fftshift(fft(u2))/N);
c  = uf(N+1:2*N);

%figure(1);  clf 
plot(x,usave(:,1),'r');  hold on;  plot(x,usave(:,1),'k.');
axis([0 2*pi -1.5 1.5]);

%  PDE solve using ODE solve

Xode = 'X29';  method = 'ode15s';
options = odeset('reltol',1e-3,'abstol',1e-6,'refine','1','stats','on');

disp(' ')
disp([' spectral diffusion: ' Xode])
disp(' ')

for tt = 0:dt:T-dt

   IVs = c;
   [t y] = feval(method,Xode,[tt tt+dt],IVs,options);
   disp(' ')

%  reconstruct u 

   c  = y(end,:)';
   uf = [0 ; -flipud(c(2:end)) ; c];
   u2 = ifft(fftshift(-i*uf))*N;
   u  = real(u2(1:N));

% interval plotting
   usave = [usave [u ; u(1)]];
   plot(x,[u ; u(1)]);

   pause(0.1)
end

plot(xg,u,'r');  hold on;  plot(xg,u,'k.');

xlabel('\bf x-axis');  ylabel('\bf u-axis')
title('\bf diffusion')