function skeleton4 ( xi, nu )

%*****************************************************************************80
%
%% skeleton4 solves a version of the 1D Allen-Cahn equation.
%
%  Discussion:
%
%    This version of the equation has the form:
%
%      du/dt = nu d2udx2 - u(u^2-1)/(2*xi)
%
%    Here nu is a positive diffusion coefficient, and xi is a parameter 
%    called the interface thickness.
%
%    The variable u is allowed the range -1 <= u <= +1.  The long-term
%    behavior of the solution is to separate into areas with u=+1 and
%    u=-1, with an intermediate interface whose width is related to xi.
%
%    The problem is solved on a 1D interval.
%
%    The boundary conditions on both ends are du/dn = 0.
%    This is handled by setting u(1)=u(2), and u(nx)=u(nx-1).
%
%    The parameter xi is described as an interface thickness.
%
%    For this problem, if nu=1, you may observe three humps appear and then:
%      xi = 0.01;   all humps stay;
%      xi = 0.015;  one hump goes away;
%      xi = 0.0155; two humps slowly disappear at different times;
%      xi = 0.02;   two humps disappear fast.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Author:
%
%    John Burkardt, Catalin Trenchea
%
%  Modified:
%
%    28 June 2020
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real XI, the interface thickness.  Default is 0.015.
%
%    real NU, the diffusion coefficient.  Default is 1.0.
%    For stability in this code, NU <= 1 is recommended.
%
  if ( nargin < 1 )
    xi = 0.015;
  end
%
%  Set the diffusion coefficient.
%
  if ( nargin < 2 )
    nu = 1.0;
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'skeleton4 \n' );
  fprintf ( 1, '  Allen-Cahn 1D reaction-diffusion equation.\n' );
  fprintf ( 1, '  Interface thickness xi = %g\n', xi );
  fprintf ( 1, '  Diffusion coefficient nu = %g\n', nu );
%
%  Set the X values.
%
  nx = 101;
  xmin = 0.0;
  xmax = 2.0;
  x = linspace ( xmin, xmax, nx );
  dx = ( xmax - xmin ) / ( nx - 1 );
%
%  Set the T values.
%
  nt = 1501;
  tmin = 0.0;
  tmax = 1.0;
  tmax = 0.15;
  t = linspace ( tmin, tmax, nt );
  dt = ( tmax - tmin ) / ( nt - 1 );
  
  fprintf ( 1, '  Use %d nodes, with dx=%g in space interval [%g,%g]\n', ...
    nx, dx, xmin, xmax );
  fprintf ( 1, '  Use %d nodes, with dt=%g in time interval [%g,%g]\n', ...
    nt, dt, tmin, tmax );
%
%  Time loop.
%
  for j = 1 : nt

    if ( j == 1 )
      u = initial_condition ( x );
    else
      uxx = laplacian_interval ( u, dx );
      dudt = ?
      u = u + dudt * dt;
    end
%
%  Apply the boundary conditions.
%
    u(1) = u(2);
    u(nx) = u(nx-1);
%
%  Save the profile at step 50.
%
    if ( j == 50 )
      u50 = u;
    end
%
%  After step 50, only plot every 10 steps.
%
    if ( j <= 50 | mod ( j, 10 ) == 1 )

      figure ( 1 );
      clf ( );
      hold ( 'on' );
      plot ( x, u, 'b-' );
      plot ( [x(1),x(nx)], [0.0,0.0], 'k-', 'linewidth', 2 ); 
      if ( 50 < j )
        plot ( x, u50, 'ro' );
        legend ( 'u(t)', 'u(step=50)' );
      end
      hold ( 'off' );
      axis ( [ xmin, xmax, -1.5, 1.5 ] );
      grid ( 'on' );
      label = sprintf ( '\\xi=%g, \\nu=%g, step=%d, t=%.3f', xi, nu, j, t(j) );
      title ( label );

      if ( j <= 50 )
        pause ( 0.2 );
      end

    end

  end
%
%  Save last plot.
%
  filename = 'exercise4.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'exercise4\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end
function u0 = initial_condition ( x )

%*****************************************************************************80
%
%% initial_condition evaluates the initial condition.
%
%  Modified:
%
%    28 June 2020
%
  xmin = min ( x );
  xmax = max ( x );

  x0 = linspace ( xmin, xmax, 84 );
  v0 = [ ...
    -1.40, -1.40, -1.44, -1.42, -1.42, ...
    -1.44, -1.43, -1.43, -1.42, -1.42, ...
    -1.40, -1.40, -1.25, -1.20, -1.17, ...
    -1.15, -1.10, -1.08, -1.00, -0.95, ...
    -0.90, -0.85, -0.88, -0.60,  0.00, ...
     0.50, -0.92, -0.25,  0.80, -0.70, ...
     0.58,  0.75,  0.58, -0.63, -0.59, ...
     0.69, -0.72,  0.70, -0.59, -0.50, ... 
     0.70, -0.79, -0.87, -0.88,  0.00, ...
     0.72, -0.80,  0.81,  0.00, -0.89, ...
     0.00,  0.70,  0.55,  0.68, -0.49, ...
     0.79,  0.00, -0.10, -0.80, -0.78, ...
    -0.83,  0.69, -0.80,  0.68,  0.50, ...
     0.70,  0.59,  1.00,  1.08,  1.10, ...
     1.15,  1.17,  1.20,  1.25,  1.30, ...
     1.30,  1.25,  1.24,  1.30,  1.31, ...
     1.30,  1.32,  1.30,  1.30 ];
%
%  Use a cubic spline to interpolate the data to the points x.
%
  u0 = csapi ( x0, v0, x );  

  return
end
function l = laplacian_interval ( u, dx )

%*****************************************************************************80
%
%% laplacian_interval approximates the laplacian on an interval.
%
%  Discussion:
%
%    The domain is represented by nx equally spaced points.
%
%    The laplacian is computed for the interior points, but a value
%    of 0 is returned for the first and last points.
%
%    It is up to the user to decide whether to reset l(1) and l(nx)
%    to nonzero values.
%
%  Modified:
%
%    28 June 2020
%
%  Input:
%
%    real u(nx): the values of a function on an equally spaced grid of points.
%
%    real dx: the spacing between any two consecutive nodes.
%
%  Output:
%
%    real l(nx): the estimate for the laplacian.
%
  nx = length ( u );
  l = zeros ( size ( u ) );
%
%  You can write this as a one-line function:
%
  l(2:nx-1) = ?;
%
%  or as a loop:
%
  for i = 2 : nx - 1
    l(i) = ?
  end
  return
end

