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

%*****************************************************************************80
%
%% kdv_etdrk4() solves the Korteweg-de Vries equation using the ETD RK4 method.
%
%  Discussion:
%
%    The system being solved is:
%
%    ut + u ux + uxxx = 0 on -pi < x < pi
%
%    or, equivalently,
%
%    ut = -1/2 d/dx(u^2) - uxxx
%
%  Modified:
%
%    13 April 2020
%
%  Author:
%
%    Lloyd Trefethen
%
%  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 X(NX): the spatial grid.
%
%    real TT(NT): the time values.
%
%    real UU(NX,NT): solution values in space and time.
%
  c8_i = sqrt ( -1.0 );
%
%  Spatial grid: -pi < x < pi, omitting last value.
%
  x = linspace ( -pi, +pi, nx + 1 );
  x = x';
  x = x(1:nx);
%
%  Set up two-soliton initial data.
%
  c1 = 25.0;
  c2 = 16.0;

  u = 3.0 * c1^2 * sech ( 0.5 * ( c1 * ( x + 2.0 ) ) ) .^ 2 ...
    + 3.0 * c2^2 * sech ( 0.5 * ( c2 * ( x + 1.0 ) ) ) .^ 2; 
  v = fft ( u );
%
%  Time step
%
  dt = 0.4 / nx^2;
%
%  Wave numbers;
%
  k = [ 0:nx/2-1, 0, -nx/2+1:-1 ]';
%
%  Fourier multipliers.
%  Note that E and E2 are defined differently in kdv_ift!
%
  L = c8_i * k .^3;
  E = exp ( dt * L );
  E2 = exp ( dt * L / 2.0 );
%
%  Roots of unity.
%
  m = 64;
  r = exp ( 2.0 * c8_i * pi * ( (1:m) - 0.5 ) / m );
%
%  Note row vector + column vector operation!
%
  LR = dt * L(:,ones(m,1)) + r(ones(nx,1),:);
%
%  Estimate Q, f1, f2, f3 by contour integral.
%
  Q  = dt * ( mean ( ...
    ( exp ( LR / 2.0 ) - 1.0 ) ./ LR, 2 ) );
  f1 = dt * ( mean ( ...
     ( - 4.0 - LR + exp ( LR ) .* ( 4.0 - 3.0 * LR + LR.^2 ) ) ./ LR.^3, 2 ) );
  f2 = dt * ( mean ( ...
    ( 2.0 + LR + exp ( LR ) .* ( - 2.0 + LR ) ) ./ LR.^3, 2 ) );
  f3 = dt * ( mean ( ...
    ( - 4.0 - 3.0 * LR - LR.^2 + exp ( LR ) .* ( 4.0 - LR ) ) ./ LR.^3, 2 ) );
%
%  Time stepping.
%
  uu = u;
  tt = 0.0;
  tmax = 0.006;
  nmax = round ( tmax / dt );
  nplt = floor ( ( tmax / 25.0 ) / 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

