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