function [uu,tt] = kursiv(x,u,dt,tmax)

% Function adaptation of kursiv.m from
% https://doi.org/10.1137/S1064827502410633
% Solves Kuramoto-Sivashinsky equation 
%     u_t = -u*u_x -u_xx -u_xxxx
% with periodic BCs by ETDRK4 scheme.
%
% Example:
%     N = 1024; 
%     x = 32*pi*(1:N)'/N; 
%     u = cos(x) + 0.1*cos(x/16).*(1+2*sin(x/16));
%     dt = 1/4;
%     tmax = 150;
%     [uu,tt] = kursiv(x,u,dt,tmax);

%% Preliminaries
N = numel(x);
assert(numel(u)==N,'mismatch between x and u');
v = fft(u); 
% Precompute various ETDRK4 scalar quantities: 
k = [0:N/2-1,0,-N/2+1:-1]'/16; % wave numbers 
L = k.^2-k.^4; % Fourier multipliers 
E = exp(dt*L); E2 = exp(dt*L/2); 
M = 16; % no. of points for complex means 
r = exp(1i*pi*((1:M)-.5)/M); % roots of unity 
LR = dt*L(:,ones(M,1)) + r(ones(N,1),:); 
Q = dt*real(mean((exp(LR/2)-1)./LR,2)); 
f1 = dt*real(mean((-4-LR+exp(LR).*(4-3*LR+LR.^2))./LR.^3,2)); 
f2 = dt*real(mean((2+LR+exp(LR).*(-2+LR))./LR.^3,2)); 
f3 = dt*real(mean((-4-3*LR-LR.^2+exp(LR).*(4-LR))./LR.^3,2)); 

%% Main time-stepping loop 
uu = u;
tt = 0; 
nmax = round(tmax/dt); 
g = -0.5i*k; 
for n = 1:nmax 
    t = n*dt; 
    Nv = g.*fft(real(ifft(v)).^2); 
    a = E2.*v+Q.*Nv; 
    Na = g.*fft(real(ifft(a)).^2); 
    b = E2.*v+Q.*Na; 
    Nb = g.*fft(real(ifft(b)).^2); 
    c = E2.*a+Q.*(2*Nb-Nv); 
    Nc = g.*fft(real(ifft(c)).^2); 
    v = E.*v+Nv.*f1+2*(Na+Nb).*f2+Nc.*f3; 
    u = real(ifft(v)); 
    uu = [uu,u]; 
    tt = [tt,t]; 
end