void rk(real a,real b,real(*f)(real,real*,real*)
	,real *y0,int n,real *x,real **y,int k
	,real h,real eps,real acc,real *gerr,real *w1,real *w2){
int i; real ymax,tol,err;

k=1; x[1]=a; *gerr=0; for(i=1,i<n,i++) y[1][i]=y0[i];    ! initialize
while(1>0){                                           ! main loop
if(x[k]+h>b) h=b-x[k];
rkstep(n,x[k],y[k],y[k+1],f,h,&err,w1,w2);
ymax=absmax(y[k+1],n); tol=ymax*eps+acc*sqrt(h/(b-a));    ! tolerance
if ( err<sqrt(2.)*tol ) {                           ! accept step
	k=k+1; x[k]=x[k-1]+h; *gerr=*gerr+err**2           ! advance x
	if(x(k)>=b) break                              ! 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)
ymax=fabs(y[k+1][1]);
for(i=2,i<n,i++) if(fabs(y[k+1][i])>ymax)ymax=fabs(y[k+1][i]
r=abs(y(1)); do i=2,n; r=max(r,abs(y(i))); end do
absmax=r; return; end