function [ xx, tt, uplt ] = allencahn_euler ( nx )

%*****************************************************************************80
%
%% allencahn_euler() solves the Allen-Cahn equation using the Euler 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.
%
%    This code is based on p34.m.
%
%  Modified:
%
%    11 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 NX: the number of nodes.
%
%  Output:
%
%    real XX(NX+1): the spatial grid
%
%    real TT(NT): the time values.
%
%    real UPLT(NT,NX): the solution values.
%

%
%  Get D, the Chebyshev differentiation matrix, and x, the Chebyshev points.
%
  [ D, x ] = cheb ( nx );
%
%  Squaring gives us the Chebyshev second derivative matrix.
%
  D2 = D^2;
%
%  Zero out the first and last rows.
%
  D2(1,:) = zeros(1,nx+1);
  D2(nx+1,:) = zeros(1,nx+1);
%
%  Initial condition.
%
  u = 0.53 * x + 0.47 * sin ( - 1.5 * pi * x );
%
%  Set the time step.
%
  eps = 0.01;
  dt = min ( [ 0.01, 50*nx^(-4)/eps ] );
  t = 0.0; 
%
%  Solve PDE by Euler formula.
%
  tmax = 100.0;
  tplot = 2;
  nplots = round ( tmax / tplot );
  plotgap = round ( tplot / dt );
  dt = tplot / plotgap;

  xx = -1:.025:1;
  uu = polyval ( polyfit ( x, u, nx ), xx );
  uplt = [ uu; zeros(nplots,length(xx)) ]; 
  tt = t;

  for i = 1 : nplots
    for n = 1 : plotgap
      t = t + dt; 
      u = u + dt * ( eps * D2 * ( u - x ) + u - u.^3 );
    end
    uu = polyval ( polyfit ( x, u, nx ), xx );
    uplt(i+1,:) = uu;
    tt = [ tt; t ];
  end

  return
end

