parameter(nmax=5,kmax=1999,pi=3.1415927)
dimension x(kmax),y0(nmax),y(nmax,kmax),work(3+6*nmax),iwork(5)
external f; common ncalls

a=0; b=4.5*pi; n=2; ncalls=0
y0(1)=0; y0(2)=1; h=0.02; eps=0.001; acc=0.001

call pc(a,b,f,n,nmax,y0,x,y,k,kmax,h,eps,acc,err,work)

print *,'#m=0, S=4'
do i=1,k ; print *,x(i),y(1,i) ; end do
print *,'#m=4, S=0'
do i=1,k ; print *,x(i),sin(x(i)) ; end do

open(1,file='out')
exact=sin(b)
write(1,*) 'required acc, eps=:',acc,eps
write(1,*) '         rk result:',y(1,k)
write(1,*) '      exact result:',exact
write(1,*) ' estimated acc,eps:',err,abs( err/y(1,k) )
write(1,*) '    actual acc,eps:',exact-y(1,k),(exact-y(1,k))/exact
write(1,*) '       evaluations:',ncalls,' rejected:',ncalls-2*(k-1)

ncalls=0; iflag=+1; !eps=0.0001; acc=0.0001
call  rkf45(f,n,y0,a,b,eps,acc,iflag,work,iwork)

write(1,*) '      rkf45 result:',y0(1)
write(1,*) '       evaluations:',ncalls
write(1,*) '             iflag:',iflag,' (should be 2)'
close(1)
end

subroutine f(x,y,dy); dimension y(*),dy(*); common ncalls
	ncalls=ncalls+1
	dy(1)=y(2); dy(2)=-y(1)
return; end