root_scalar_by_type Subroutine

    subroutine root_scalar_by_type(method,fun,ax,bx,xzero,fzero,iflag,&
                                   ftol,rtol,atol,maxiter,fax,fbx,&
                                   bisect_on_failure)

    implicit none

    type(root_method),intent(in)  :: method   !! the method to use
    procedure(func2)              :: fun      !! user function to find the root of
    real(wp),intent(in)           :: ax       !! left endpoint of initial interval
    real(wp),intent(in)           :: bx       !! right endpoint of initial interval
    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 (`-1`=error, `0`=root found, `-999`=invalid method)
    real(wp),intent(in),optional  :: ftol     !! absolute tolerance for `f=0`
    real(wp),intent(in),optional  :: rtol     !! relative tol for x
    real(wp),intent(in),optional  :: atol     !! absolute tol for x
    integer,intent(in),optional   :: maxiter  !! maximum number of iterations
    real(wp),intent(in),optional  :: fax      !! if `f(ax)` is already known, it can be input here
    real(wp),intent(in),optional  :: fbx      !! if `f(bx)` is already known, it can be input here
    logical,intent(in),optional   :: bisect_on_failure  !! if true, then if the specified method fails,
                                                        !! it will be retried using the bisection method.
                                                        !! (default is False). Note that this can use up
                                                        !! to `maxiter` additional function evaluations.

    class(root_solver),allocatable :: s

    select case (method%id)

    case(root_method_brent%id);                allocate(brent_solver                :: s)
    case(root_method_bisection%id);            allocate(bisection_solver            :: s)
    case(root_method_regula_falsi%id);         allocate(regula_falsi_solver         :: s)
    case(root_method_illinois%id);             allocate(illinois_solver             :: s)
    case(root_method_anderson_bjorck%id);      allocate(anderson_bjorck_solver      :: s)
    case(root_method_ridders%id);              allocate(ridders_solver              :: s)
    case(root_method_pegasus%id);              allocate(pegasus_solver              :: s)
    case(root_method_bdqrf%id);                allocate(bdqrf_solver                :: s)
    case(root_method_muller%id);               allocate(muller_solver               :: s)
    case(root_method_brenth%id);               allocate(brenth_solver               :: s)
    case(root_method_brentq%id);               allocate(brentq_solver               :: s)
    case(root_method_chandrupatla%id);         allocate(chandrupatla_solver         :: s)
    case(root_method_toms748%id);              allocate(toms748_solver              :: s)
    case(root_method_zhang%id);                allocate(zhang_solver                :: s)
    case(root_method_anderson_bjorck_king%id); allocate(anderson_bjorck_king_solver :: s)
    case(root_method_blendtf%id);              allocate(blendtf_solver              :: s)
    case(root_method_barycentric%id);          allocate(barycentric_solver          :: s)
    case(root_method_itp%id);                  allocate(itp_solver                  :: s)

    case default
        iflag = -999    ! invalid method
        return
    end select

    call s%initialize(func_wrapper,ftol,rtol,atol,maxiter)
    call s%solve(ax,bx,xzero,fzero,iflag,fax,fbx,bisect_on_failure)

    contains

        function func_wrapper(me,x) result(f)
            implicit none
            class(root_solver),intent(inout) :: me
            real(wp),intent(in) :: x
            real(wp) :: f
            f = fun(x)
        end function func_wrapper

    end subroutine root_scalar_by_type