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

a=-10; b=10; n=1; ncalls=0
y0(1)=0; ; h=0.5; eps=0.001; acc=0.001

call rk(a,b,f2,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),atan(x(i))-atan(-10.) ; end do

open(1,file='out')
exact=atan(10.)-atan(-10.)
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)=1./(1.+x*x)
return; end

subroutine f2(n,x,y,dy); dimension y(n),dy(n); common ncalls
	ncalls=ncalls+1
	dy(1)=1./(1.+x*x)
return; end