function [ x, tt, uu ] = kursiv_pde_etdrk4 ( nx )

%*****************************************************************************80
%
%% kursiv_pde_etdrk4() solves the Kuramoto-Sivashinsky PDE using the ETD RK4 method.
%
%  Discussion:
%
%    The system being solved is:
%
%    ut = - u * ux - uxx - uxxx, periodic BC on 0 < x < 32 pi
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    12 April 2020
%
%  Author:
%
%    Original MATLAB version by Lloyd Trefethen.
%    This version by John Burkardt.
%
%  Reference:
%
%    Stephen Cox, Paul Matthews,
%    Exponential time differencing for stiff systems,
%    Journal of Computational Physics,
%    Volume 176, pages 430-455, 2002.
%
%    Aly-Khan Kassam, Lloyd Trefethen,
%    Fourth-order time-stepping for stiff ODE's,
%    SIAM Journal on Scientific Computing,
%    Volume 26, Number 4, pages 1214-1233, 2005.
%
%    Lloyd Trefethen,
%    Spectral methods in MATLAB,
%    SIAM, 2000,
%    LC: QA377.T65
%    ISBN: 978-0-898714-65-4
%
%  Input:
%
%    integer NX: the number of nodes.
%
%  Output:
%
%    real XX(NX+1): the spatial grid
%
%    real TT(NT): the time values.
%
%    real UU(NX,NT):
%
  c8_i = sqrt ( -1.0 );
%
%  Spatial grid: 0 < x < 32 * pi
%
  x = 32.0 * pi * ( 1 : nx )' / nx;
%
%  Initial condition.
%
  u = cos ( x / 16.0 ) .* ( 1.0 + sin ( x / 16.0 ) );
  v = fft ( u );
%
%  Precompute various ETD RK4 scalar quantities.
%
  dt = 1.0 / 4.0;
  k = [ 0:nx/2-1, 0.0, -nx/2+1:-1 ]' / 16.0;
  L = k.^2 - k.^4;
  E = exp ( dt * L );
  E2 = exp ( dt * L / 2.0 );
  M = 16;
%
%  Roots of unity.
%
  r = exp ( c8_i * pi * ( (1:M) - 0.5 ) / M );
%
%  Note row vector + column vector operation creates matrix LR!
%
  LR = dt * L(:,ones(M,1)) + r(ones(nx,1),:);

  Q  = dt * real ( mean ( ...
    ( exp ( LR / 2.0 ) - 1.0 ) ./ LR, 2 ) );
  f1 = dt * real ( mean ( ...
     ( - 4.0 - LR + exp ( LR ) .* ( 4.0 - 3.0 * LR + LR.^2 )) ./ LR.^3, 2 ) );
  f2 = dt * real ( mean ( ...
    ( 2.0 + LR + exp ( LR ) .* ( - 2.0 + LR ) ) ./ LR.^3, 2 ) );
  f3 = dt * real ( mean ( ...
    ( - 4.0 - 3.0 * LR - LR.^2 + exp ( LR ) .* ( 4.0 - LR ) ) ./ LR.^3, 2 ) );
%
%  Time stepping.
%
  uu = u;
  tt = 0.0;
  tmax = 150.0;
  nmax = round ( tmax / dt );
  nplt = floor ( ( tmax / 100 ) / dt );
  g = - 0.5 * c8_i * k;

  for i = 1 : nmax

    t = i * 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.0 * Nb - Nv );
    Nc = g .* fft ( real ( ifft ( c ) ) .^2 );

    v = E .* v + Nv .* f1 + 2.0 * ( Na + Nb ) .* f2 + Nc .* f3;

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

  end

  uu = uu';

  return
end

