function y = reaction_exact ( t )

%*****************************************************************************80
%
%% reaction_exact() returns the exact solution of reaction_ode().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    23 January 2021
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Willem Hundsdorfer, Jan Verwer,
%    Numerical solution of time-dependent advection-diffusion-reaction equations,
%    Springer, 2003
%    ISBN: 978-3-662-09017-6
%
%  Input:
%
%    real T: the current time.
%
%  Output:
%
%    real Y(3): the exact solution.
%
  [ k, t0, y0, tstop ] = reaction_parameters ( );

  a0 = y0(1);
  b0 = y0(2);
  c0 = y0(3);
  d0 = a0 - b0;

  if ( d0 == 0 )
    qt = k * t;
  else 
    qt = ( 1.0 - exp ( - k * t * d0 ) ) / d0;
  end

  at = a0 ./ ( 1.0 + b0 * qt );
  bt = at - d0;
  ct = c0 + ( a0 - at );

  y = [ at, bt, ct ];

  return
end
