function [ x, tt, uu ] = burgers_pde_etdrk4 ( nx, nt, vis )

%*****************************************************************************80
%
%% burgers_pde_etdrk4() solves the Burgers equation using the ETD RK4 method.
%
%  Discussion:
%
%    The system being solved is:
%
%      ut = -1/2 d/dx(u^2) + vis * uxx
%
%    over
%
%      -pi < x < +pi
%
%    with initial condition
%
%      u = exp ( - 10.0 * ( sin ( 0.5 * x ) ).^2 );
%
%    and periodic boundary conditions.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    20 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.
%    A typical value is NX = 500.
%
%    integer NT: the number of time steps to be plotted.
%
%    real VIS: the viscosity.  
%    "Interesting" values are VIS = 0.03 and 0.
%
%  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.
%
  x = linspace ( -pi, pi, nx + 1 );
  x = x';
  x = x(1:nx);
%
%  Initial condition.
%
  u = exp ( - 10.0 * ( sin ( 0.5 * x ) ).^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.
%
  L = c8_i * vis * k .^ 2;
  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 creates a matrix LR!
%
  LR = dt * L(:,ones(m,1)) + r(ones(nx,1),:);
%
%  Estimate Q, f1, f2, f3 by contour integral.
%
  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.
%
  tmax = 1.0;
  nmax = round ( tmax / dt );
  jstep = floor ( nmax / ( nt - 1 ) );

  uu = u;
  tt = 0.0;

  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, jstep ) == 0 )
      u = real ( ifft ( v ) );
      uu = [ uu, u ];
      tt = [ tt; t ];
    end

  end

  uu = uu';

  return
end

