c
c  @(#)rkqc.f	1.1
c
	SUBROUTINE RKQC(Y,DYDX,N,X,HTRY,EPS,YSCAL,HDID,HNEXT,DERIVS)
	PARAMETER (NMAX=10,PGROW=-0.20,PSHRNK=-0.25,FCOR=1./15.,
     *		ONE=1.,SAFETY=0.9,ERRCON=6.E-4)
	EXTERNAL DERIVS
	DIMENSION Y(N),DYDX(N),YSCAL(N),YTEMP(NMAX),YSAV(NMAX),DYSAV(NMAX)
	XSAV=X
	DO I=1,N
		YSAV(I)=Y(I)
		DYSAV(I)=DYDX(I)
	ENDDO
	H=HTRY
1	HH=0.5*H
	CALL RK4(YSAV,DYSAV,N,XSAV,HH,YTEMP,DERIVS)
	X=XSAV+HH
	CALL DERIVS(X,YTEMP,DYDX)
	CALL RK4(YTEMP,DYDX,N,X,HH,Y,DERIVS)
	X=XSAV+H
	IF(X.EQ.XSAV)PAUSE 'Stepsize not significant in RKQC.'
	CALL RK4(YSAV,DYSAV,N,XSAV,H,YTEMP,DERIVS)
	ERRMAX=0.
	DO I=1,N
		YTEMP(I)=Y(I)-YTEMP(I)
		ERRMAX=MAX(ERRMAX,ABS(YTEMP(I)/YSCAL(I)))
	ENDDO
	ERRMAX=ERRMAX/EPS
	IF(ERRMAX.GT.ONE) THEN
		H=SAFETY*H*(ERRMAX**PSHRNK)
		GOTO 1
	ELSE
		HDID=H
		IF(ERRMAX.GT.ERRCON)THEN
			HNEXT=SAFETY*H*(ERRMAX**PGROW)
		ELSE
			HNEXT=4.*H
		ENDIF
	ENDIF
	DO I=1,N
		Y(I)=Y(I)+YTEMP(I)*FCOR
	ENDDO
	RETURN
	END
