#include <math.h>

#define PGROW -0.20
#define PSHRNK -0.25
#define FCOR 0.06666666		/* 1.0/15.0 */
#define SAFETY 0.9
#define ERRCON 6.0e-4

int rkqc(y,dydx,x,htry,eps,yscal,hdid,hnext,derivs)
     double y[],dydx[],*x,htry,eps,yscal[],*hdid,*hnext;
     void (*derivs)();	/* ANSI: void (*derivs)(double,double *,double *); */
{
  double xsav,hh,h,temp,errmax;
  double dysav[6],ysav[10],ytemp[10];
  extern void rk4();
  extern double Cinc();

  xsav=(*x);

  ysav[0]=y[0];
  ysav[1]=y[1];
  ysav[2]=y[2];
  ysav[3]=y[3];
  ysav[4]=y[4];
  ysav[5]=y[5];
  dysav[0]=dydx[0];
  dysav[1]=dydx[1];
  dysav[2]=dydx[2];
  dysav[3]=dydx[3];
  dysav[4]=dydx[4];
  dysav[5]=dydx[5];

  h=htry;
  ytemp[6]=ysav[6]=y[6];  ytemp[7]=ysav[7]=y[7];  
  ytemp[8]=ysav[8]=y[8];  ytemp[9]=ysav[9]=y[9];

  for (;;) {
    hh=0.5*h;
    rk4(ysav,dysav,xsav,hh,ytemp,derivs);
    *x=xsav+Cinc(hh);
    (*derivs)(dydx,(int) *x,ytemp);
    rk4(ytemp,dydx,*x,hh,y,derivs);
    *x=xsav+Cinc(h);
    if (*x == xsav) 
      return(0); /* Convergence  Failure. */
    rk4(ysav,dysav,xsav,h,ytemp,derivs);
    errmax=0.0;

    ytemp[0] -= y[0];
    if (errmax < (temp=fabs(ytemp[0]/yscal[0]))) errmax=temp;
    ytemp[1] -= y[1];
    if (errmax < (temp=fabs(ytemp[1]/yscal[1]))) errmax=temp;
    ytemp[2] -= y[2];
    if (errmax < (temp=fabs(ytemp[2]/yscal[2]))) errmax=temp;
    ytemp[3] -= y[3];
    if (errmax < (temp=fabs(ytemp[3]/yscal[3]))) errmax=temp;
    ytemp[4] -= y[4];
    if (errmax < (temp=fabs(ytemp[4]/yscal[4]))) errmax=temp;
    ytemp[5] -= y[5];
    if (errmax < (temp=fabs(ytemp[5]/yscal[5]))) errmax=temp;

    if ((errmax /= eps) <= 1.0) {
      *hdid=h;
      *hnext=(errmax > ERRCON ? SAFETY*h*pow(errmax, PGROW) : 4.0*h);
      break;
    }
    h=SAFETY*h*exp(PSHRNK*log(errmax));
  }
  
  y[0] += ytemp[0]*FCOR;
  y[1] += ytemp[1]*FCOR;
  y[2] += ytemp[2]*FCOR;
  y[3] += ytemp[3]*FCOR;
  y[4] += ytemp[4]*FCOR;
  y[5] += ytemp[5]*FCOR;
  return(1); /* success */
}

#undef PGROW
#undef PSHRNK
#undef FCOR
#undef SAFETY
#undef ERRCON