double rkqs(double[] y,double[] dydt,double t,double htry,double[] yscal){
  int i;
  double errmax, h, htemp, tnew;
  h=htry;
  for (; ; ) {
    rkck(y,dydt,t,h);
    errmax=0.0;
    for (i=0; i < nvar; i++)     errmax=Math.max(errmax,Math.abs(Yerr[i] / yscal[i]));
    errmax/=eps;
    if (errmax <= 1.0)     break;
    if (Double.isNaN(errmax)) {
      throw new RuntimeException("errmax NaN");
    }
 else {
      htemp=SAFETY * h * Math.pow(errmax,PSHRNK);
      h=(h >= 0.0 ? Math.max(htemp,0.1 * h) : Math.min(htemp,0.1 * h));
    }
    tnew=t + h;
    if (tnew == t) {
      throw new RuntimeException("stepsize underflow in rkqs");
    }
  }
  if (errmax > ERRCON)   hnext=SAFETY * h * Math.pow(errmax,PGROW);
 else   hnext=5.0 * h;
  t+=(hdid=h);
  for (i=0; i < nvar; i++)   y[i]=Ytemp[i];
  return t;
}

double rkqs ( double [ ] y , double [ ] dydt , double t , double htry , double [ ] yscal ) { int i ; double errmax , h , htemp , tnew ; h = htry ; for ( ; ; ) { rkck ( y , dydt , t , h ) ; errmax = 0.0 ; for ( i = 0 ; i < nvar ; i ++ ) errmax = Math . max ( errmax , Math . abs ( Yerr [ i ] / yscal [ i ] ) ) ; errmax /= eps ; if ( errmax <= 1.0 ) break ; if ( Double . isNaN ( errmax ) ) { throw new RuntimeException ( "errmax NaN" ) ; } else { htemp = SAFETY * h * Math . pow ( errmax , PSHRNK ) ; h = ( h >= 0.0 ? Math . max ( htemp , 0.1 * h ) : Math . min ( htemp , 0.1 * h ) ) ; } tnew = t + h ; if ( tnew == t ) { throw new RuntimeException ( "stepsize underflow in rkqs" ) ; } } if ( errmax > ERRCON ) hnext = SAFETY * h * Math . pow ( errmax , PGROW ) ; else hnext = 5.0 * h ; t += ( hdid = h ) ; for ( i = 0 ; i < nvar ; i ++ ) y [ i ] = Ytemp [ i ] ; return t ; }
take an adaptive step
