function [t,u,e] = ssrkc4(neqn,t0,u0,h)
      [ut0]=derv(neqn,t0,u0);
      k1=h*ut0;
      u=u0+0.5*k1;
      t=t0+0.5*h;
      [ut]=derv(neqn,t,u);
      k2=h*ut;
      u=u0+0.5*k2;
      t=t0+0.5*h;
      [ut]=derv(neqn,t,u);
      k3=h*ut;
      u=u0+k3;
      t=t0+h;
      [ut]=derv(neqn,t,u);
      k4=h*ut;
      u=u0+k2;
      sum2=u0+(1.0/6.0)*(k1+2.0*k2+2.0*k3+k4);
      t=t0+h;
      e=sum2-u;
      u=u+e;
