function rkf(x0,y0,h,xf,epsilon,epsb,sigma)
%
%  driver for the Runge-Kutta Fehlberg method
%
j=1;
Y(j)=y0;
x(j)=x0;
nfun=0;
while (x(j) < xf)
   [Yhat,Yn]=comprkf(x(j),Y(j),h);
   nfun=nfun+6;
   while (abs(Yn-Yhat) > epsilon)
      h=h/2;
      [Yhat,Yn]=comprkf(x(j),Y(j),h);
      nfun=nfun+6;
   end
   if (abs(Yn-Yhat) <= epsb) 
      h=sigma*h;
   end
   Y(j+1)=Yn;
   x(j+1)=x(j)+h;
   j=j+1;
end
plot(x,Y)
nfun

