function y = flame_exact ( t )

%*****************************************************************************80
%
%% flame_exact evaluates the exact solution of the flame ODE.
%
%  Discussion:
%
%    1 equation.
%
%    Moler attributes this problem to Lawrence Shampine.
%
%    The equation describes the radius of a ball of flame that
%    begins, at time 0, at DELTA.
%
%      Y(0) = DELTA
%
%    The rate of fuel consumption is proportional to the volume, and
%    the rate of fuel intake is proportional to the area of the ball.
%    We take the constant of proportionality to be 1.
%
%      Y' = Y^2 - Y^3
%
%    The data is normalized so that Y = 1 is the equilibrium solution.
%
%    The computation is to be made from T = 0 to T = 2/DELTA.
%
%    For values of DELTA close to 1, such as 0.01, the equation is
%    not stiff.  But for DELTA = 0.0001, the equation can become
%    stiff as the solution approaches the equilibrium solution Y = 1,
%    and computed solutions may be wildly inaccurate or cautious
%    solvers may take very small timesteps.
%
%    The exact solution involves the Lambert W function, defined by
%
%      W(z) * exp ( W(z) ) = z
%
%    and if we set
%
%      A = ( 1 / DELTA - 1 )
%
%    then
%
%      Y(T) = 1 / ( W(A*exp(A-T)) + 1 )
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    10 November 2020
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Cleve Moler,
%    Cleve's Corner: Stiff Differential Equations,
%    MATLAB News and Notes,
%    May 2003, pages 12-13.
%
%  Input:
%
%    real T(:): the times.
%
%  Output:
%
%    real Y(:), the exact solution values.
%
  [ t0, y0 ] = flame_parameters ( );

  a = ( 1.0 - y0 ) / y0;
  y = 1.0 ./ ( lambertw ( a .* exp ( a - ( t - t0 ) ) ) + 1.0 );

  return
end
