      function [u]=rkf45a(neqn,t0,tf,u0,npts)
      h=(tf-t0)/npts;
      for i=1:npts
         [t,u,e]=ssrkf45(neqn,t0,u0,h);
         u0=u;
         t0=t;
      end
