subroutine rkstep(n,x,y0,y1,f,h,err,z0,z12)
dimension y0(n),y1(n),z0(n),z12(n)

call f(x,y0,z0)
do i=1,n
	y1(i)=y0(i)+.5*h*z0(i)
enddo
call f(x+h/2,y1,z12)
do i=1,n
	y1(i)=y0(i)+h*z12(i)
enddo
err=0
do i=1,n
	dy=h*(z0(i)-z12(i))/6
	err=max(err,abs(dy))
end do
return; end