function [ f1, f2, f3, f4, f5, s ] = r8_fehl ( f, neqn, y, t, h, yp )

%*****************************************************************************80
%
%% R8_FEHL takes one Fehlberg fourth-fifth order step.
%
%  Discussion:
%
%    This routine integrates a system of NEQN first order ordinary differential
%    equations of the form
%      dY(i)/dT = F(T,Y(1:NEQN))
%    where the initial values Y and the initial derivatives
%    YP are specified at the starting point T.
%
%    The routine advances the solution over the fixed step H and returns
%    the fifth order (sixth order accurate locally) solution
%    approximation at T+H in array S.
%
%    The formulas have been grouped to control loss of significance.
%    The routine should be called with an H not smaller than 13 units of
%    roundoff in T so that the various independent arguments can be
%    distinguished.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 August 2010
%
%  Author:
%
%    Original FORTRAN77 version by Herman Watts, Lawrence Shampine.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Erwin Fehlberg,
%    Low-order Classical Runge-Kutta Formulas with Stepsize Control,
%    NASA Technical Report R-315, 1969.
%
%    Lawrence Shampine, Herman Watts, S Davenport,
%    Solving Non-stiff Ordinary Differential Equations - The State of the Art,
%    SIAM Review,
%    Volume 18, pages 376-411, 1976.
%
%  Parameters:
%
%    Input, function F, a user-supplied function to evaluate the
%    derivatives dYdT, of the form:
%
%      function yp = f ( t, y )
%
%    Input, integer NEQN, the number of equations to be integrated.
%
%    Input, real Y(NEQN), the current value of the
%    dependent variable.
%
%    Input, real T, the current value of the independent
%    variable.
%
%    Input, real H, the step size to take.
%
%    Input, real YP(NEQN), the current value of the
%    derivative of the dependent variable.
%
%    Output, real F1(NEQN), F2(NEQN), F3(NEQN), F4(NEQN),
%    F5(NEQN), derivative values needed for the computation.
%
%    Output, real S(NEQN), the estimate of the solution at T+H.
%
  ch = h / 4.0;

  f5(1:neqn) = y(1:neqn) + ch * yp(1:neqn);

  f1 = f ( t + ch, f5 );

  ch = 3.0 * h / 32.0;

  f5(1:neqn) = y(1:neqn) + ch * ( yp(1:neqn) + 3.0 * f1(1:neqn) );

  f2 = f ( t + 3.0 * h / 8.0, f5 );

  ch = h / 2197.0;

  f5(1:neqn) = y(1:neqn) + ch * ...
  ( 1932.0 * yp(1:neqn) ...
  + ( 7296.0 * f2(1:neqn) - 7200.0 * f1(1:neqn) ) ...
  );

  f3 = f ( t + 12.0 * h / 13.0, f5 );

  ch = h / 4104.0;

  f5(1:neqn) = y(1:neqn) + ch * ...
  ( ...
    ( 8341.0 * yp(1:neqn) - 845.0 * f3(1:neqn) ) ...
  + ( 29440.0 * f2(1:neqn) - 32832.0 * f1(1:neqn) ) ...
  );

  f4 = f ( t + h, f5 );

  ch = h / 20520.0;

  f1(1:neqn) = y(1:neqn) + ch * ...
  ( ...
    ( -6080.0 * yp(1:neqn) ...
    + ( 9295.0 * f3(1:neqn) - 5643.0 * f4(1:neqn) ) ...
    ) ...
  + ( 41040.0 * f1(1:neqn) - 28352.0 * f2(1:neqn) ) ...
  );

  f5 = f ( t + h / 2.0, f1 );
%
%  Ready to compute the approximate solution at T+H.
%
  ch = h / 7618050.0;

  s(1:neqn) = y(1:neqn) + ch * ...
  ( ...
    ( 902880.0 * yp(1:neqn) ...
    + ( 3855735.0 * f3(1:neqn) - 1371249.0 * f4(1:neqn) ) ) ...
  + ( 3953664.0 * f2(1:neqn) + 277020.0 * f5(1:neqn) ) ...
  );

  return
end