subroutine quindibisect2(c,b,d,n,m1,m2,eps1,relfeh,eps2,z,x) 
        real (kind=8) :: eps1,eps2,relfeh,xmin,xmax,h,x1,xu,xo
        real (kind=8),dimension(n):: c,b,d,x
        real (kind=8) ,dimension(:):: wu(m1:m2)
        integer :: i,n,m1,m2,z,a,k
        integer maxit
        parameter (maxit=1000)
        
        d(1)=0
        d(2)=0
        b(1)=0
        xmin = c(n)-dabs(b(n))-dabs(d(n))
        xmax = c(n) + dabs(b(n)) + dabs(d(n))
        h=dabs(b(n-1))+dabs(d(n-1))+dabs(b(n))
        if ((c(n-1)+h).gt.xmax) xmax=c(n-1)+h
        if ((c(n-1)-h).lt.xmin) xmin=c(n-1)-h
        do i=n-2,1,-1
            h= dabs(b(i))+dabs(d(i))+dabs(b(i+1))+dabs(d(i+2))
            if ((c(i)+h).gt.xmax) xmax=c(i)+h
            if ((c(i)-h).lt.xmin) xmin=c(i)-h
        enddo
        if ((xmin+xmax).gt.0) then 
           eps2=relfeh*xmax
        else 
           eps2=-relfeh*xmin
        endif
        if (eps1.lt.0) eps1=eps2
        eps2= 0.5*eps1+7*eps2
        xo=xmax
        do i=m1,m2 
          x(i)=xmax
          wu(i)=xmin
        enddo
        z=0
        do k=m2,m1,-1
          xu=xmin
!comment loop for the kth value
          do i=k,m1,-1
            if (xu.lt.wu(i)) then 
              xu=wu(i) 
              go to 10
            endif
          enddo 
       10 if (xo.gt.x(k)) xo=x(k)
          x1=(xu+xo)/2
          do while  (xo-xu.gt.2*relfeh*(dabs(xu)+dabs(xo)+eps1))
             if (z.gt.maxit) exit
             z=z+1
             call quindet2(c,b,d,n,x1,a)
             a=n-a
             if (a.lt.k) then
                if (a.lt.m1) then 
                   wu(m1)=x1
                   xu=x1
                else
                  wu(a+1)=x1
                  xu=x1
                  if (x(a).gt.x1) x(a)=x1
                endif
             else
                xo=x1
             endif
             x1=(xu+xo)/2
          enddo
          x(k)=(xo+xu)/2
        enddo ! k
       endsubroutine