function [t,u,e] = ssrkf45(neqn,t0,u0,h)
         [ut0]=derv(neqn,t0,u0);
         k1=h*ut0;
         u=u0+0.25*k1;
         t=t0+0.25*h;
         [ut]=derv(neqn,t,u);
         k2=h*ut;
         u=u0+(3.0/32.0)*k1...
             +(9.0/32.0)*k2;
         t=t0+(3.0/8.0)*h;
         [ut]=derv(neqn,t,u);
         k3=h*ut;
         u=u0+(1932.0/2197.0)*k1...
             -(7200.0/2197.0)*k2...
             +(7296.0/2197.0)*k3;
         t=t0+(12.00/13.0)*h;
         [ut]=derv(neqn,t,u);
         k4=h*ut;
         u=u0+( 439.0/ 216.0)*k1...
             -(   8.0       )*k2...
             +(3680.0/ 513.0)*k3...
             -( 845.0/4104.0)*k4;
         t=t0+h;
         [ut]=derv(neqn,t,u);
         k5=h*ut;
         u=u0-(   8.0/  27.0)*k1...
             +(   2.0       )*k2...
             -(3544.0/2565.0)*k3...
             +(1859.0/4104.0)*k4...
             -(  11.0/  40.0)*k5;
         t=t0+0.5*h;
         [ut]=derv(neqn,t,u);
         k6=h*ut;
         u=u0+(    25.0/ 216.0)*k1...
             +(  1408.0/2565.0)*k3...
             +(  2197.0/4104.0)*k4...
             -(     1.0/   5.0)*k5;
      sum2=u0+(   16.0/  135.0)*k1...
             +( 6656.0/12825.0)*k3...
             +(28561.0/56430.0)*k4...
             -(    9.0/   50.0)*k5...
             +(    2.0/   55.0)*k6;
         t=t0+h;
         e=sum2-u;
         u=u+e;
