function [ xx, tt, uu ] = allencahn_etdrk4 ( nx )

%*****************************************************************************80
%
%% allencahn_etdrk4() solves the Allen-Cahn equation using the ETD RK4 method.
%
%  Discussion:
%
%    The system being solved is:
%
%    ut = 0.01 * uxx + u - u^3 for -1 < x < +1
%    u(-1) = -1
%    u(+1) = +1
%
%    The computation is based on Chebyshev points, so the linear term
%    is nondiagonal.
%
%  Modified:
%
%    11 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 (minus 1).
%
%  Output:
%
%    real XX(NX+1): the spatial grid.
%
%    real TT(NT): the time values.
%
%    real UU(NT,NX+1): the solution values.
%
  c8_i = sqrt ( -1.0 );
%
%  Set parameter.
%
  alpha = 0.01;
%
%  Get the Chebyshev differentiation matrix and grid of NX+1 points.
%
  [ D, x ] = cheb ( nx );
%
%  Drop first and last points, leaving NX-1 points.
%
  x = x(2:nx);
%
%  Initial condition.
%  Let w = u - x, so boundary conditions become homogeneous.
%
  w = 0.53 * x + 0.47 * sin ( - 1.5 * pi * x ) - x;
  u = [ 1.0; w+x; -1.0 ];
%
%  Time step.
%
  dt = 1.0 / 4.0;
%
%  Precompute various ETD RK4 matrix quantities.
%
  M = 32;
  r = 15.0 * exp ( c8_i * pi * ( (1:M) - 0.5 ) / M );
  L = D^2;
  L = alpha * L(2:nx,2:nx);
  A = dt * L;
  E = expm ( A );
  E2 = expm ( A / 2.0 );
%
%  Estimate Q, f1, f2, f3 with a contour integral.
%
  I = eye ( nx - 1 );
  Z = zeros ( nx - 1 );

  Q  = Z;
  f1 = Z;
  f2 = Z;
  f3 = Z;

  for j = 1 : M
    z = r(j);
    zIA = inv ( z * I - A );
    Q  = Q +  dt * zIA * ( exp ( z / 2.0 ) - 1.0 );
    f1 = f1 + dt * zIA * ( - 4.0 - z + exp ( z ) * ( 4.0 - 3.0 * z + z^2 ) ) / z^2;
    f2 = f2 + dt * zIA * ( 2.0 + z + exp ( z ) * ( z - 2.0 ) ) / z^2;
    f3 = f3 + dt * zIA * ( - 4.0 - 3.0 * z - z^2 + exp ( z ) * ( 4.0 - z ) ) / z^2;
  end

  Q  = real ( Q  / M );
  f1 = real ( f1 / M );
  f2 = real ( f2 / M );
  f3 = real ( f3 / M );
%
%  Time stepping.
%
  uu = u;
  tt = 0.0;
  tmax = 70.0;
  nmax = round ( tmax / dt );
  nplt = floor ( ( tmax / 70 ) / dt );

  for i = 1 : nmax

    t = i * dt;

    Nu = ( w + x ) - ( w + x ) .^ 3;
    a = E2 * w + Q * Nu;
    Na = a + x - ( a + x ) .^ 3;
    b = E2 * w + Q * Na;
    Nb = b + x - ( b + x ) .^ 3;
    c = E2 * a + Q * ( 2.0 * Nb - Nu );
    Nc = c + x - ( c + x ) .^ 3;

    w = E * w + f1 * Nu + 2.0 * f2 * ( Na + Nb ) + f3 * Nc;

    if ( mod ( i, nplt ) == 0 )
      u = [ 1.0; w+x; -1.0 ];
      uu = [ uu, u ];
      tt = [ tt, t ];
    end

  end

  uu = uu';
  xx = [ 1.0; x; -1.0 ];

  return
end

