function [U,x,t] = kuramotoShivashinskyIntegrate(u,Lx,dt,Nt)

% Crank-Nicolson/Adams-Bashforth (CNAB2) integration of the 1D
% Kuramoto-Shivashinsky equation 
%     $u_t = -uu_x - u_{xx} - u_{xxxx}$
% on [0,L] with periodic boundary conditions, using Fourier modes for
% spatial discretization and Crank-Nicolson/Adams-Bashforth (CNAB2)
% timestepping for temporal discretization. Adapted/elaborated from 
%     https://online.kitp.ucsb.edu/online/transturb17/gibson/html/5-kuramoto-sivashinksy.html
%
% Example:
%     Lx = 128;
%     Nx = 1024;
%     dt = 1/16;
%     Nt = 1600;
%     %
%     x = Lx*(0:(Nx-1))/Nx;
%     u = cos(x) + 0.1*cos(x/16).*(1+2*sin(x/16));
%     %
%     [U,x,t] = kuramotoShivashinskyIntegrate(u,Lx,dt,Nt);
%     figure; pcolor(x,t,real(U)); shading interp;
%     xlabel('x'); ylabel('t');
%
% For background, see
%     https://en.wikipedia.org/wiki/Crank%E2%80%93Nicolson_method#The_method
%     https://en.wikipedia.org/wiki/Linear_multistep_method#Two-step_Adams%E2%80%93Bashforth
% CN is used for the linear term; AB2 (i.e, two-step AB) is used for the
% nonlinear term. Their synthesis for the KS equation is described in 
%     https://scicomp.stackexchange.com/questions/37336/
%     https://online.kitp.ucsb.edu/online/transturb17/gibson/html/5-kuramoto-sivashinksy.html
% and this code is adapted from the latter source.
% 
% The KS equation can be written
%     $u_t = Lu + Nu$
% where $Lu := -u_{xx} - u_{xxxx}$ and $Nu := -uu_x$ are respectively linear
% and nonlinear terms. 
%
% The CNAB2 scheme is 
%     $[u(x,t+dt) - u(x,t)]/dt = 
%         0.5 L[u(x,t+dt) + u(x,t)] + 1.5 Nu(x,t) - 0.5 Nu(x,t-dt)$
% which can be rewritten as 
%     $(I - 0.5 dt L)u(x,t+dt) = (I + 0.5 dt L)u(x,t)
%         1.5 dt Nu(x,t) - 0.5 dt Nu(x,t-dt)$.
% In Fourier space, $L$ is diagonal: writing $A := I + 0.5 dt L)$ and $B :=
% (I - 0.5 dt L)^{-1}$, the above becomes 
%     $u(x,t+dt) = B[Au(x,t)+1.5 dt Nu(x,t) - 0.5 dt Nu(x,t-dt)]$
%
% Last modified 20211019 by Steve Huntsman

%% Preliminaries
Nx = length(u);                 % number of gridpoints
kx = [0:Nx/2-1,0,-Nx/2+1:-1];	% integer wavenumbers: exp(2*pi*i*kx*x/L)
alpha = 2*pi*kx/Lx;             % real wavenumbers:    exp(i*alpha*x)
D = 1i*alpha;                   % D = d/dx operator in Fourier space
L = -D.^2-D.^4;

% Convenient diagonal linear operators (see top comments)
A = ones(1,Nx)+0.5*dt*L;
B = (ones(1,Nx)-0.5*dt*L).^(-1);

% Spectral representation of nonlinear term
NuCur = -0.5*D.*fft(u.*u);      % current time step

%% Initialization
x = (0:(Nx-1))*Lx/Nx;
t = (0:(Nt-1))*dt;
U = nan(Nt,Nx);
U(1,:) = u;                     % initial condition
u = fft(u);                     % convert to spectral representation

%% Main loop
for n = 1:Nt
    % Nonlinear term: shift in time and update
	NuPre = NuCur;
	NuCur  = -0.5*D.*fft(real(ifft(u)).^2);
    
    % CNAB2 (see top comments)
	u = B.*(A.*u+1.5*dt*NuCur-0.5*dt*NuPre);
    U(n,:) = fft(u);
end