function [ x, tdata, udata ] = kdv_ift ( nx )

%*****************************************************************************80
%
%% kdv_ift() solves the Korteweg-deVries equation using the IFT method.
%
%  Discussion:
%
%    The system being solved is:
%
%      ut + u ux + uxxx = 0 on -pi < x < pi
%
%    by FFT with integrating factor v = exp(-ik^3t)*uhat.
%
%    This code is related to program p27.m in the Trefethen reference.
%
%  Modified:
%
%    12 April 2020
%
%  Author:
%
%    Lloyd Trefethen
%
%  Reference:
%
%    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 N: the number of nodes.
%
%  Output:
%
%    real X(NX): the spatial grid
%
%    real TDATA(NT): the time values.
%
%    real UDATA(NX,NT): solution values at space and time grid points.
%
  c8_i = sqrt ( -1.0 );
%
%  Set up grid.
%
  dt = 0.4 / nx^2;
  x = ( 2.0 * pi / nx ) * ( - nx / 2 : nx / 2 - 1 )';
%
%  Set up two-soliton initial data.
%
  A = 25.0;
  B = 16.0;

  u = 3.0 * A^2 * sech ( 0.5 * ( A * ( x + 2.0 ) ) ) .^ 2 ...
    + 3.0 * B^2 * sech ( 0.5 * ( B * ( x + 1.0 ) ) ) .^ 2; 
  v = fft ( u );
%
%  Define the Fourier multipliers.
%  Note that E and E2 are different in kdv_etdrk4!
%
  k = [ 0:nx/2-1, 0, -nx/2+1:-1 ]';
  ik3 = c8_i * k .^ 3;
  E = exp ( dt * ik3 / 2.0 );
  E2 = E .^ 2;
%
%  Solve PDE.
%
  tmax = 0.006;
  nplt = floor ( ( tmax / 25.0 ) / dt );
  nmax = round ( tmax / dt );
  udata = u;
  tdata = 0.0;

  g = - 0.5 * c8_i * dt * k;

  for i = 1 : nmax

    t = i * dt;

    a = g .* fft ( real ( ifft (         v             ) ) .^ 2 );
    b = g .* fft ( real ( ifft ( E  .* ( v + a / 2.0 ) ) ) .^ 2 );
    c = g .* fft ( real ( ifft ( E  .*   v + b / 2.0   ) ) .^ 2 );
    d = g .* fft ( real ( ifft ( E2 .*   v + E .* c    ) ) .^ 2 );

    v = E2 .* v + ( E2 .* a + 2.0 * E .* ( b + c ) + d ) / 6.0;

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

  end

  udata = udata';

  return
end

