c
c  @(#)rk4.f	1.1
c
	SUBROUTINE RK4(Y,DYDX,N,X,H,YOUT,DERIVS)
	PARAMETER (NMAX=10)
	DIMENSION Y(N),DYDX(N),YOUT(N),YT(NMAX),DYT(NMAX),DYM(NMAX)
	HH=H*0.5
	H6=H/6.
	XH=X+HH
	DO I=1,N
		YT(I)=Y(I)+HH*DYDX(I)
	ENDDO
	CALL DERIVS(XH,YT,DYT)
	DO I=1,N
		YT(I)=Y(I)+HH*DYT(I)
	ENDDO
	CALL DERIVS(XH,YT,DYM)
	DO I=1,N
		YT(I)=Y(I)+H*DYM(I)
		DYM(I)=DYT(I)+DYM(I)
	ENDDO
	CALL DERIVS(X+H,YT,DYT)
	DO I=1,N
		YOUT(I)=Y(I)+H6*(DYDX(I)+DYT(I)+2.*DYM(I))
	ENDDO
	RETURN
	END
