c*******************************************************
      SUBROUTINE BISECT (xa,xb,nmax,xr,errfun)
c*******************************************************
c     roots of a nonlinear function FF  -  Bisection method
c
      IMPLICIT REAL*8 (a-h, o-z)
c
c-- check of starting data
      fa=FF(xa)
      fb=FF(xb)
      IF (dabs(fa).lt.errfun) THEN  !  You're lucky, this is the solution
         xr=xa
         errfun=dabs(fa)
         return
      ELSEIF (dabs(fb).lt.errfun) THEN  !  You're lucky, this is the solution
         xr=xb
         errfun=dabs(fb)
         return
      ELSEIF (fa*fb.gt.0.d0) THEN  !  Impossible to go futher
         write(*,*) fa,fb
         STOP  'error data (BISECT):  f(xa) * f(xb) > 0'
      ENDIF
c
c-- Convergence history recorded in file CONV
        open(1,file='conv')
        rewind(1)
        write(1,*)xa,fa
        write(1,*)xb,fb
c
c___ start iterating _______________
c
      DO niter=1,nmax
         xr=0.5d0*(xa+xb)
         fr=FF(xr)
         write(1,*)xr,fr
         write(*,*) niter,xr,fr    !   screen message
c-- Is XR a root?
         IF (dabs(fr) .lt. errfun) THEN     !  Yes, it is ....
            errfun=dabs(fr)
            goto 10
         ELSE                               !  No, it isn't
c-- define the new interval
         IF (fa*fr .LT. 0.d0) THEN  !  continue with the LEFT sub-interval
            xb=xr
            fb=fr
         ELSE                       !  continue with the RIGHT sub-interval
            xa=xr
            fa=fr
         ENDIF
         ENDIF
      ENDDO
        stop   'Max. Number of Iterations reached WITHOUT CONVERGENCE '
 10   close(1)
      RETURN
      END

      include 'testfun.for'