ModAB Subroutine

    subroutine ModAB(me,ax,bx,fax,fbx,xzero,fzero,iflag)

    implicit none

    class(modab_solver),intent(inout) :: me
    real(wp),intent(in)    :: ax      !! left endpoint of initial interval
    real(wp),intent(in)    :: bx      !! right endpoint of initial interval
    real(wp),intent(in)    :: fax     !! `f(ax)`
    real(wp),intent(in)    :: fbx     !! `f(ax)`
    real(wp),intent(out)   :: xzero   !! abscissa approximating a zero of `f` in the interval `ax`,`bx`
    real(wp),intent(out)   :: fzero   !! value of `f` at the root (`f(xzero)`)
    integer,intent(out)    :: iflag   !! status flag (`0`=root found, `-2`=max iterations reached)

    real(wp) :: x1,x2,x3,y1,y2,y3,ym,m,r,k
    integer :: i  !! iteration counter
    logical :: root_found, bis
    integer :: side !! for tracking the side
    real(wp) :: threshold !! threshold to fall back to bisection if AB fails to shrink the interval enough

    iflag = 0
    side = 0
    x1 = ax; y1 = fax
    x2 = bx; y2 = fbx
    bis = .true.
    threshold = x2-x1
    do i = 1, me%maxiter
        if (bis) then
            x3 = 0.5_wp*(x1+x2)
            y3 = me%f(x3)
            if (me%solution(x3,y3,xzero,fzero)) return
            ym = 0.5_wp*(y1+y2)
            r  = 1.0_wp - abs(ym/(y2-y1))  ! symmetry factor
            k  = r*r                       ! deviation factor
            if (abs(ym-y3) < k*(abs(ym) + abs(y3))) then ! function is close enough to linear
                bis = .false.
                threshold = 16.0_wp*(x2-x1)   ! safety factor of 4 bisection iters = 2^4
            end if
        else
            x3 = (x1*y2-y1*x2)/(y2-y1)
            ! Clamp secant point to the interval to handle floating-point errors
            if (x3 <= x1) then
                x3 = x1
                y3 = y1  ! clamped: reuse known y1
            else if (x3 >= x2) then
                x3 = x2
                y3 = y2  ! clamped: reuse known y2
            else
                y3 = me%f(x3)  ! not clamped: evaluate y3
                if (me%solution(x3,y3,xzero,fzero)) return
            end if
            threshold = 0.5_wp * threshold
        end if

        ! convergence check:
        root_found = me%converged(x1,x2)
        if (root_found .or. i==me%maxiter) then
            xzero = x3
            fzero = y3
            if (.not. root_found) iflag = -2  ! max iterations reached
            exit
        end if
        select case (side)
        case(1)
            m = 1.0_wp - y3/y1
            if (m<=0.0_wp) then
                y2 = y2 * 0.5_wp
            else
                y2 = y2 * m
            end if
        case(2)
            m = 1.0_wp - y3/y2
            if (m<=0.0_wp) then
                y1 = y1 * 0.5_wp
            else
                y1 = y1 * m
            end if
        end select

        if (sign(1.0_wp,y1) == sign(1.0_wp,y3)) then
            if (.not. bis) side = 1
            x1 = x3
            y1 = y3
        else
            if (.not. bis) side = 2
            x2 = x3
            y2 = y3
        end if
        if (x2-x1 > threshold) then ! if Anderson-Bjork is not shrinking the interval fast enough
            bis  = .true. ! reset to bisection.
            side = 0
        end if

    end do

    end subroutine ModAB