barycentric Subroutine

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

        implicit none

        class(barycentric_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)

        integer :: i !! iteration counter
        real(wp) :: k !! a real number in the interval (0, 0.5)
        real(wp) :: x0,x1,x2,x3,x4,x5,x6,&
                    f0,f1,f2,f3,f4,f5,f6,&
                    w0,w1,w2,w3,w22,w32

        real(wp),parameter :: bisection_tol = 1.0e-17_wp  !! if `f0` or `f1` are below this value,
                                                          !! the method switches to bisection.
                                                          !! (should be an input or computed from `epsilon`?)
        real(wp),parameter :: k0 =  1.0_wp / 6.0_wp !! initial value of `k` (should be an input?)

        iflag = 0
        x0 = ax; f0 = fax
        x1 = bx; f1 = fbx
        k = k0

        do i = 1, me%maxiter

            if (i>1) then
                f0 = me%f(x0)
                if (me%solution(x0,f0,xzero,fzero)) return
                f1 = me%f(x1)
                if (me%solution(x1,f1,xzero,fzero)) return
            end if

            if (me%converged(x0, x1) .or. i == me%maxiter) then
                if (i == me%maxiter) iflag = -2 ! max iterations reached
                exit
            end if

            if (abs(f0)<bisection_tol .or. abs(f1)<bisection_tol) then
                ! if f is too small, switch to bisection
                k = 0.5_wp
            end if

            x2 = x0 + k*(x1-x0)   ! first point for interpolation
            f2 = me%f(x2)
            if (me%solution(x2,f2,xzero,fzero)) return

            if (f0*f2 < 0.0_wp) then
                x1 = x2
            else if (f0 == f2) then
                x0 = x2
            else
                x3 = x1 + k*(x0-x1) ! second point for interpolation
                f3 = me%f(x3)
                if (me%solution(x3,f3,xzero,fzero)) return

                if (f3*f1 < 0.0_wp) then
                    x0 = x3
                else if (f3 == f1) then
                    x0 = x2        ! note: typo in python version here?
                    x1 = x3
                else

                    w0 = 1.0_wp/(f0*(f0-f2)*(f0-f3))   ! first quadratic interpolation
                    w2 = 1.0_wp/(f2*(f2-f0)*(f2-f3))
                    w3 = 1.0_wp/(f3*(f3-f0)*(f3-f2))
                    x4 = (w0*x0+w2*x2+w3*x3)/(w0+w3+w2)
                    f4 = me%f(x4)
                    if (me%solution(x4,f4,xzero,fzero)) return

                    w1  = 1.0_wp/(f1*(f1-f2)*(f1-f3))  ! second quadratic interpolation
                    w22 = 1.0_wp/(f2*(f2-f1)*(f2-f3))
                    w32 = 1.0_wp/(f3*(f3-f1)*(f3-f2))
                    x5  = (w1*x1+w22*x2+w32*x3)/(w1+w32+w22)
                    f5  = me%f(x5)
                    if (me%solution(x5,f5,xzero,fzero)) return

                    ! determine if x3,x4 are in the interval [x2,x3]
                    if (x4<x2 .or. x4>x3) then
                        if (x5<x2 .or. x5>x3) then
                            x0 = x2
                            x1 = x3
                        else
                            if (f5*f2<0.0_wp) then
                                x0 = x2
                                x1 = x5
                            else
                                x0 = x5
                                x1 = x3
                            end if
                        end if
                    else
                        if (x5<x2 .or. x5>x3) then
                            if (f4*f2<0.0_wp) then
                                x0 = x2
                                x1 = x4
                            else
                                x0 = x4
                                x1 = x3
                            end if
                        else

                            if (f5*f4>0.0_wp) then
                                if (f4*f2<0.0_wp) then
                                    x0 = x2
                                    x1 = min(x4,x5)
                                else
                                    x0 = max(x4,x5)
                                    x1 = x3
                                end if
                            else
                                ! cubic interpolation if f4 and f5 have opposite signs
                                w0 = w0/(f0-f1)
                                w1 = w1/(f1-f0)
                                w2 = w2/(f2-f1)
                                w3 = w3/(f3-f1)
                                x6 = (w0*x0+w1*x1+w2*x2+w3*x3)/(w0+w1+w3+w2)
                                f6 = me%f(x6)
                                if (me%solution(x6,f6,xzero,fzero)) return

                                ! verify that x6 is in the interval [x2,x3]
                                if (x6<=x2 .or. x6>=x3) then
                                    x0 = min(x4,x5)
                                    x1 = max(x4,x5)
                                else
                                    if (f6*f4<0.0_wp) then
                                        x0 = min(x4,x6)
                                        x1 = max(x4,x6)
                                    else
                                        x0 = min(x5,x6)
                                        x1 = max(x5,x6)
                                    end if
                                end if
                            end if
                        end if
                    end if
                end if
            end if

        end do

        call choose_best(x0,x1,f0,f1,xzero,fzero)

    end subroutine barycentric