% p32.m - Solve Kuramoto-Sivashinsky(KS) eq. 
%           u_t + uu_x + u_xx + u_xxxx = 0 on [-pi L,pi L] by
%         FFT with integrating factor.

clf; clear;

function [x tdata udata] = KSsolu(uex, L, N, tmax, dt, nplot) 
  x = (2*pi*L/N)*(0:N-1)';
  u = feval(uex, x);    v = fft(u); 
  k = [0:N/2-1 0 -N/2+1:-1]'/L;
  
  nplt = floor((tmax/nplot)/dt); nmax = round(tmax/dt);
  udata = u; tdata = 0; 
  for n = 1:nmax
    t = n*dt; g = -.5i*dt*k;
    E = exp(-dt*(k.^4-k.^2)/2); E2 = E.^2;
    a = g.*fft(real( ifft(     v    ) ).^2);
    b = g.*fft(real( ifft(E.*(v+a/2)) ).^2);     % 4th-order
    c = g.*fft(real( ifft(E.*v + b/2) ).^2);     % Runge-Kutta
    d = g.*fft(real( ifft(E2.*v+E.*c) ).^2);
    v = E2.*v + (E2.*a + 2*E.*(b+c) + d)/6;
    if mod(n, nplt) == 0
       u = real(ifft(v)); 
       udata = [udata u]; tdata = [tdata t];
    end
  end
end

L=16; u =inline('cos(x/L).*(1+sin(x/L))', 'x'); 

tmax=300; dt=0.05; nplot=300;  N=128;
[x, tdata, udata] = KSsolu(u, L, N, tmax, dt, nplot);  

[tt, xx]=meshgrid(tdata,x);  imagesc(tt, xx, udata); xlabel t, ylabel x




