subroutine pc(a,b,f,n,nmax,y0,x,y,k,kmax,h,eps,acc,gerr,w)
dimension x(kmax),y(nmax,kmax),y0(nmax),w(5*nmax)
external f

k=1; x(1)=a; gerr=0; do i=1,n; y(i,1)=y0(i); end do    ! initialization
10 continue;                                           ! main loop
if(x(k)+h>b) h=b-x(k)
if(k==1) then                                          ! first step
	call rkstep(n,x(k),y(1,k),y(1,k+1),f,h,err,w(1),w(n+1))
else
	call pcstep(x(k-1),x(k),h,y(1,k-1),y(1,k),y(1,k+1),n,f,err,&
&	w(1),w(n+1),w(2*n+1),w(3*n+1))
end if
ymax=absmax(y(1,k+1),n); tol=(ymax*eps+acc)*sqrt(h/(b-a))
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.25*0.9)
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