% m42_kuramotosivashinsky.m - Kuramoto-Sivashinsky equation
%
%         u_t = -u_xx - u_xxxx - (u^2/2)_x
%
% on [-16pi,16pi] with periodic boundary conditions.
% Long waves grow because of -u_xx;
% short waves decay because of -u_xxxx; 
% the nonlinear term transfers energy from long to short.

% Grid, initial data, and plotting setup:
  npts = 400;
  h = 32*pi/npts;
  x = -16*pi + (1:npts)'*h;
  u = cos(x/16).*(1+sin(x/16));
  hold off, shg
  plt = plot(x,u,'linewidth',2);
  axis([-16*pi 16*pi -4 4]), grid
  k = 0.2; tmax = input('tmax? (e.g. 100) ');
  disp('type <return> to see solution'), pause

% First sparse matrix for the linear, implicit term:
  L = length(x);
  a = 1 + 6*k/h^4 - 2*k/h^2; b = -4*k/h^4 + k/h^2; c = k/h^4;
  main = a*sparse(ones(L,1));
  off  = b*sparse(ones(L-1,1));
  off2 = c*sparse(ones(L-2,1));
  B = diag(main) + diag(off,1) + diag(off,-1) ...
                 + diag(off2,2) + diag(off2,-2);
  B(end,1) = b; B(1,end) = b;
  B(end-1,1) = c; B(end,2) = c; B(1,end-1) = c; B(2,end) = c;

% Second sparse matrix for the nonlinear, explicit term:
  d = k/(2*h);
  off  = d*sparse(ones(L-1,1));
  A = diag(off,1) - diag(off,-1);
  A(end,1) = d; A(1,end) = -d;
  
% Time-stepping:
  t = 0;                  % Key feature of this discretization:
  while t < tmax-1e-6     %   implicit for linear terms (so k can be big)
    u = B\(u - A*u.^2/2); %   explicit for nonlinear term (so no nonlin eqs)
    t = t+k; set(plt,'ydata',u)
    title(['t = ' num2str(t)],'fontsize',20)
    drawnow
  end