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