c***  4th order runge-kutta, more or less straight out of 
c***  numerical recipes

	subroutine rk4(y,dydx,par,n,x,h,yout,derivs)

	implicit none

	integer i, n
	real y(n),dydx(n),yout(n)
	real yt(n), par(n), dyt(n)
	real dym(n)
	real h, hh, x, xh, h6
	external derivs

	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,par,n)
	do i=1,n
	 yt(i)=y(i)+hh*dyt(i)
	enddo
	call derivs(xh,yt,dym,par,n)
	do i=1,n
	 yt(i)=y(i)+h*dym(i)
	 dym(i)=dyt(i)+dym(i)
	enddo
	call derivs(x+h,yt,dyt,par,n)
	do i=1,n
	 yout(i)=y(i)+h6*(dydx(i)+dyt(i)+2.*dym(i))
	enddo

	return
	end