function [u]=rkc4b(neqn,t0,tf,u0,npts,abserr,relerr)
global HINTEG
h=(tf-t0)/2;
hmin=(tf-t0)/npts;
t=t0;
while t <= tf*0.999 
   if h < hmin
      h = hmin;
      u0=u;
      t0=t;
   end 
   if t+h > tf
      h=tf-t;
   end                                        
   [t,u,e]=ssrkc4(neqn,t0,u0,h);
   nfin1=1;
   for i=1:neqn
      if abs(e(i)) > abs(u(i))*relerr+abserr
         nfin1=0;
         h=h/2.0;
         break;
      end
   end
   if nfin1 == 1
      u0=u;
      t0=t;
      for i=1:neqn
         if abs(e(i)) > (abs(u(i))*relerr+abserr)/16.
            nfin1=0;
            break;
         end
      end
      if nfin1 == 1
         h=h*2.0;
         if t+h > tf
            h=tf-t;
         end 
      end
   end
end