function [tout, yout] = sir54(FunFcn, t0, tend, y0, tol);
% y'= FunFcn(x,y)
% y(t0) = y0
% interval [t0,tend]
%
% SIR example with beta=0.25, gamma=0.1, popul=1e7, begin with 1 infection
% >> fcn=@(x,y)[-0.25*y(1)*y(2)/1e7;0.25*y(1)*y(2)/1e7-0.1*y(2);gamma*y(2)];
% >> [tout,yout] = sir54(fcn,0,200,[9999999 1 0]',0.001);


a=[0,0,0,0,0,0,0;
   0.224991857145594237,0,0,0,0,0,0;
   0.0891533610320671103,0.238960211116888553,0,0,0,0,0;
   1.30502239863058047,-5.50177549088639896,5.14113139529316096,0,0,0,0;
   1.76234012927544337,-7.44800428470202672,6.70879999021246612, ...
   -0.0341377330357163063,0,0,0;
   1.87412911619818746,-7.91467977711770718,7.07985467092264069, ...
   -0.0244945338238275367,-0.0148094761792934300,0,0;
   0.0986971498551664256,0.00100729346874150652,0.495118366873759549, ...
   3.93767069618687179,-13.7395424087173685,10.2070489023328292,0]' ;
c=[0,0.224991857145594237,0.328113572148955663, ...
   0.944378303037342471,0.988998101750166470,1,1]';
b=[0.0986971498551664256, 0.00100729346874150652,0.495118366873759549, ...
   3.93767069618687179,-13.7395424087173685, 10.2070489023328292,0];
bb=[0.0898860401768279254,0.000893217827115906704,0.524893254061430804, ...
    2.14886680528438840, -5.7054302422968618, 3.8908909249470988,0.05];
%-------------------------------------------------------
% Initialization
t = t0;hmax = (tend - t)/5;hmin = (tend - t)/2000000;
y = y0(:);ytrue=y0;f = y*zeros(1,length(c));
tout = t;yout = y.';h=tol0.2/max(max(abs(f(:,1))),1);
f(:,7) = feval(FunFcn,t,y);ireject=0;

% The main loop
while (t < tend) & (h >= hmin)
   if t + h > tend, h = tend - t; end
   if ireject==0,f(:,1)=f(:,7);end %FSAL

   % Compute the six stages
   for j = 2:7,
       f(:,j) = feval(FunFcn, t+c(j)*h, y+h*f*a(:,j));
   end
   % Error estimation
   delta = h*norm(h*f*(b-bb)','inf');

   % Compare tolerance and error estimation and update the solution
   if delta <= tol,
      t = t + h;y = y + h*f*b';
      tout = [tout; t];
      yout = [yout; y.'];
      ireject=0;
   else
      irejct=1; % step rejected, do not use FSAL
   end

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

if (t < tend)
   disp('SINGULARITY LIKELY.');t
end