function [tout, yout,ireject] = scalar64(FunFcn, t0, tfinal, y0, tol);
%	Integrate a scalar autonomus differential equation using
%	4th and 6th order Runge-Kutta formulas.  See also ODE45
%
%	[T,Y] = scalar64('yprime', T0, Tfinal, Y0) integrates the scalar
%	autonomus ordinary differential equation described by the M-file
%	YPRIME.M over the interval T0 to Tfinal and using initial
%	condition Y0.
%
% INPUT:
% F     - String containing name of user-supplied problem description.
%         Call: yprime = fun(y) where F = 'fun'.
%         t      - Time (scalar).
%         y      - Solution (scalar).
%         yprime - Returned derivative scalar; yprime(i) = dy(i)/dt.
% t0    - Initial value of t.
% tfinal- Final value of t.
% y0    - Initial value scalar.
% tol   - The desired accuracy. (Default: tol = 1.e-6).
%
% OUTPUT:
% T  - Returned integration time points (row-vector).
% Y  - Returned solution, one solution scalar per tout-value.
% ireject - number of rejected steps
%
% The result can be displayed by: plot(tout, yout).
%
% example solve y'=-y^3/2, y(0)=1 over [0,2] (with exact solution y=1/sqrt(1+t))
% [tout, yout] = scalar64(@(y) -y^3/2, 0, 1, 1, 1e-6);
% [tout yout 1./sqrt(1+tout)]
%
%   Simos & Tsitouras, 6th order Runge-Kutta pairs for scalar autonomous
%   Appl. Comput. Math., 19 (2020) 412-421 
%-------------------------------------------------------

%---- rk6(4) pair ----
beta=[797978303/2260428898,0,0,0,0,0,0;
99311147/342166459,30757295/439108719,0,0,0,0,0;
13389565/580674501,-146017213/215544313,379248055/279391273,0,0,0,0;
99230792/216123993,280001283/177111158,-352130836/197579617,255765624/374426161,0,0,0;
2523438/281492519,-434261384/254809577,673986568/353490849,-124667480/284104839,100015407/312846799,0,0;
2519828/168294781,-96974223/380218696,129477794/219833563,87056077/279905590,22900559/152069488,22283152/117599949,0]';
alpha=[797978303/2260428898,122072623/338820474,284537668/404729185,487737511/518352431,22842351/247609048,1]';
gamma=[[ 2519828/168294781,-96974223/380218696,129477794/219833563,87056077/279905590,22900559/152069488,22283152/117599949,0]
       [ -0.003032795661621  -0.172351240781884   0.184830608078977  -0.036383060837677   0.073917105080843   0.003019384121362  -0.05]]';
%-------------------------------------------------------

ireject=0;
pow = 1/5; 
if nargin < 5, tol = 1.e-6; end

% Initialization
t = t0;
hmax = (tfinal - t)/5;
hmin = (tfinal - t)/2000000;
h = (tfinal - t)/100;
y = y0(:);
f = y*zeros(1,length(alpha)+1);
tout = t;
yout = y.';
f(:,1) = feval(FunFcn,y);
h=tol^pow/max(max(abs(f(:,1))),1);

% The main loop
   while (t < tfinal) & (h >= hmin)
      if t + h > tfinal, h = tfinal - t; end

      % Compute the slopes
      f(:,1) = feval(FunFcn,y);
      for j = 1:length(alpha),
         f(:,j+1) = feval(FunFcn, y+h*f*beta(:,j));
      end

      % Estimate the error and the acceptable error
      delta = norm(h^2*f*gamma(:,2),'inf'); 

      % Update the solution only if the error is acceptable
      if delta <= tol,
         t = t + h;
         y = y + h*f*gamma(:,1);
         tout = [tout; t];
         yout = [yout; y.'];
      else
         ireject=ireject+1; 
      end

      % Update the step size
      if delta ~= 0.0
         h = min(hmax, 0.9*h*(tol/delta)^pow);
      end
   end;

   if (t < tfinal)
      disp('SINGULARITY LIKELY.')
      t
   end