% ks_cnab2.m - solution of Kuramoto-Sivashinsky equation by CNAB2. 
% u_t = u*u_x - u_xx - u_xxxx, periodic BCs on [0,L]

% computation is based on v = fft(u), so linear term is diagonal
% compare p27.m in Trefethen, "Spectral Methods in MATLAB", SIAM 2000
% AK Kassam and LN Trefethen, July 2002

% Parameters: space, time, discretization 
Lx = 32*pi;       % length of domain
N  = 512;         % number gridpoints  
dt = 1/16;        % time step
T  = 100;         % integration time
x  = Lx*(1:N)'/N; % spatial grid

% Set initial condition
u = cos(x/16).*(1+sin(x/16)).*(1+0.01*sin(x/32));
v = fft(u);

% Precompute various CNAB2 quantities:
q = (2*pi/Lx)*[0:N/2-1 0 -N/2+1:-1]';  % wave numbers
L = q.^2 - q.^4;                       % Fourier multipliers for linear part of SH eqn
D = i*q;                               % differentiation operator
g = -0.5*D;                            % -1/2 d/dx, for -1/2 d/dx(u^2) = -u u_x

% Initialize time-stepping: 
Nv  = g.*fft(u.^2);
Nvp = Nv;

LR = (1 + (dt/2)*L);
LL = (1 - (dt/2)*L);
LL_inv = LL.^(-1);

nmax = round(T/dt); 
nplt = floor((T/100)/dt);
uu = u;
tt = 0;

for n = 1:nmax
  t = n*dt;

  Nvp = Nv;
  Nv = g.*fft(real(ifft(v)).^2);
  
  v = LL_inv.*(LR.*v + 1.5*dt*Nv - 0.5*dt*Nvp);
  %v = LL_inv.*(LR.*v);
  

  if mod(n,nplt)==0
    u = real(ifft(v));
    uu = [u, uu]; tt = [tt,t];
  end
end

%min(min(uu))
%max(max(uu))
% Plot results:
surf(tt,x,uu), shading interp, lighting phong, axis tight
view([90 90]), colormap(jet); 
set(gca,'zlim',[0 50])
%light('color',[1 1 1],'position',[-1,2,2])
%caxis([-6 2])
%material([0.30 0.60 0.60 40.00 1.00]);
%material('dull');
xlbl('t')
ylbl('x')
%set(gca, 'xtick', [])
%set(gca, 'ytick', [])
%set(gca,'xcolor',[0.8 0.8 0.8])
%set(gca,'ycolor',[0.8 0.8 0.8])
colorbar