subroutine rk(a,b,f,n,nmax,y0,x,y,k,kmax,h,eps,acc,gerr,work)
dimension x(kmax),y(nmax,kmax),y0(nmax),work(2*nmax)
external f

k=1; x(1)=a; gerr=0; do i=1,n; y(i,1)=y0(i); end do    ! initialize
10 continue; if(x(k)+h>b) h=b-x(k)                     ! main loop
call rkstep(n,x(k),y(1,k),y(1,k+1),f,h,err,work(1),work(n+1))
ymax=absmax(y(1,k+1),n)
tol=max(ymax*eps,acc*sqrt(h/(b-a)))/sqrt(2.)           ! tolerance
if ( err<sqrt(2.)*tol ) then                           ! accept step
	k=k+1; x(k)=x(k-1)+h; gerr=gerr+err**2           ! advance x
	if(x(k)>=b) goto 20                              ! done
	if(k .eq. kmax) stop 'rk2: k>kmax'
	h=h*stpest(tol,err); goto 10
else                                                   ! reject step
	h=h*stpest(tol,err); goto 10
end if
20 continue; gerr=sqrt(gerr); return; end

function stpest(tol,err)
if(err.eq.0) then
	stpest=4
else
	stpest=min(4.,abs(tol/err)**0.2)*0.95
end if
return; end

function absmax(y,n)
dimension y(n)
r=abs(y(1)); do i=2,n; r=max(r,abs(y(i))); end do
absmax=r; return; end