subroutine rkstep(n,x,y0,y1,f,h,err,z0,z12,z1)
dimension y0(n),y1(n),z0(n),z12(n),z1(n)
!
! advances solution of a system of ordinary differential equations
! dy_i/dx = f_i(x,y)
! by one step h using Runge-Kutta method
!
call f(n,x,y0,z0)
do i=1,n; y1(i)=y0(i)+h/2*z0(i); end do; call f(n,x+h/2,y1,z12)
do i=1,n; y1(i)=y0(i)+h*z0(i);   end do; call f(n,x+h,y1,z1)
do i=1,n; y1(i)=y0(i)+h*(z0(i)+4*z12(i)+z1(i))/6; end do
err=0; do i=1,n
	dy=h/2*( y1(i)-(y0(i)+h*z12(i)) )/3
	err=max(err,abs(dy))
end do
return; end