bis_secant Subroutine

    PURE SUBROUTINE bis_secant   ( error   , fnctn   , maxerr  , maxsteps,  &
                                   steps   , xcur    , xleft   , xright     )
      !! A root finder using the hybrid bisection-secant algorithm.
      !! Returns a computed value of the root within xleft and xright
      !! once error < maxerr or has run maximum number of iterations.
        IMPLICIT NONE

        interface
            pure function fnctn ( x        )
              !! The Fortran function for which the root will be
              !! found. Must take only one real argument, return a
              !! real argument.
                real(8), intent(in) ::  x
                real(8)             ::  fnctn
            end function fnctn
        end interface
        
        integer, intent(in)     ::  maxsteps
          !! An integer value for the maximum number of iterations
          !! applied to the algorithm before terminating.
        real(8), intent(in)     ::  maxerr
          !! A real value which specifies the the maximum allowable
          !! error in the computed root. The subroutine terminates
          !! once the error passes below this level.
        real(8), intent(inout)  ::  xleft
          !! A real value specifying lower bound within which to
          !! search for root.
        real(8), intent(inout)  ::  xright
          !! A real value specifying upper bound within which to
          !! search for root.
        integer, intent(out)    ::  steps
          !! Returns the number of iterations needed before error
          !! falls below maxerr (or returns maxsteps if that many
          !! iterations occur first).
        real(8), intent(out)    ::  error
          !! A real variable in which an estimate of the error in the
          !! computed root will be stored and returned.
        real(8), intent(out)    ::  xcur
          !! A real variable in which the computed value of the root
          !! will be returned.
        
        ! Other variables:
        REAL(8) ::  dfxcur,                                     &
                    fxcur,                                      &
                    fxleft,                                     &
                    fxprev1,                                    &
                    fxright,                                    &
                    xnext,                                      &
                    xprev1
    !------------------------------------------------------------------!

        ! Initialize variables
        fxleft = fnctn(xleft)
        fxright = fnctn(xright)
        dfxcur = 0.d0
        fxcur = 0.d0
        fxprev1 = 0.d0
        xnext = 0.d0
        xprev1 = 0.d0
        
        ! Check that there is a root bracketed by 'xleft' and 'xright'
        IF ( fxleft * fxright > 0.0d0 ) THEN
           error stop ('No root in specified interval.')
        END IF
        
        ! Initialize variables)
        xcur = xleft
        fxcur = fxleft
        xprev1 = xright
        fxprev1 = fxright
        
        ! Iterate until root converges or reach 'maxsteps'
        DO steps = 1, maxsteps
            ! If not 1st iteration, update variables for next iteration
            IF ( steps /= 1 ) THEN
                fxprev1 = fxcur
                fxcur = fnctn(xcur)
                IF ( fxleft * fxcur < 0.d0 ) THEN
                    xright = xcur
                    fxright = fxcur
                ELSE
                    xleft = xcur
                    fxleft = fxcur
                END IF
            END IF
            
            ! Compute a secant and use to find 'xnext'
            dfxcur = (fxcur - fxprev1) / (xcur - xprev1)
            xnext = xcur - fxcur / dfxcur
            ! If this 'xnext' value is outside of left or right bounds, 
            ! instead use a bisection to find 'xnext'
            IF ( ( xnext < xleft ) .OR. ( xnext > xright ) ) THEN
                xnext = 5.d-1 * (xleft + xright)
            END IF
            
            ! Estimate error in approximating root and update variables
            error = ABS(xnext - xcur)
            xprev1 = xcur
            xcur = xnext
            
            ! If error less than tolerance, return
            IF ( error < maxerr ) RETURN
        END DO
        
        ! If reached maximum iterations, print a warning message
        error stop ('No solution found in specified number of steps.')

    END SUBROUTINE bis_secant