1 subroutine quindibisect2(c,b,d,n,m1,m2,eps1,relfeh,eps2,z,x)
2 real (kind=8) :: eps1,eps2,relfeh,xmin,xmax,h,x1,xu,xo
3 real (kind=8),dimension(n):: c,b,d,x
4 real (kind=8) ,dimension(:):: wu(m1:m2)
5 integer :: i,n,m1,m2,z,a,k
12 xmin = c(n)-dabs(b(n))-dabs(d(n))
13 xmax = c(n) + dabs(b(n)) + dabs(d(n))
14 h=dabs(b(n-1))+dabs(d(n-1))+dabs(b(n))
15 if ((c(n-1)+h).gt.xmax) xmax=c(n-1)+h
16 if ((c(n-1)-h).lt.xmin) xmin=c(n-1)-h
18 h= dabs(b(i))+dabs(d(i))+dabs(b(i+1))+dabs(d(i+2))
19 if ((c(i)+h).gt.xmax) xmax=c(i)+h
20 if ((c(i)-h).lt.xmin) xmin=c(i)-h
22 if ((xmin+xmax).gt.0) then
27 if (eps1.lt.0) eps1=eps2
37 !comment loop for the kth value
44 10 if (xo.gt.x(k)) xo=x(k)
46 do while (xo-xu.gt.2*relfeh*(dabs(xu)+dabs(xo)+eps1))
49 call quindet2(c,b,d,n,x1,a)
58 if (x(a).gt.x1) x(a)=x1