%
%  lect01.m -- spectral diffusion w/ nonlinearity
%           -- bifurcation & instability
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

%  parameters

global n2 N D a b

choice = 1

switch choice
  case 1
    D =  1.0;  a =  3.0;   b = 0.5;
    N =   64;  L = pi/2;  dx = 2*L/N;  
    T =   32; dt =  4.0;

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  = (4/pi^2) * (xg.*(2*L-xg))';

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),'g');  hold on;  plot(x,usave(:,1),'k.');
axis([0 2*L -1.5 1.5]);

%  timestep PDE using matlab ODE solver

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

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 & nonlinearity')

data = [data ; a max(u)]
