integrate_to_event_fixed_step Subroutine

    subroutine integrate_to_event_fixed_step(me,t0,x0,h,tmax,tol,tf,xf,gf)

    implicit none

    class(rk_fixed_step_class),intent(inout) :: me
    real(wp),intent(in)               :: t0      !! initial time
    real(wp),dimension(:),intent(in)  :: x0      !! initial state
    real(wp),intent(in)               :: h       !! abs(time step)
    real(wp),intent(in)               :: tmax    !! max final time if event not located
    real(wp),intent(in)               :: tol     !! function tolerance for root finding
    real(wp),intent(out)              :: tf      !! actual final time reached
    real(wp),dimension(:),intent(out) :: xf      !! final state (at tf)
    real(wp),intent(out)              :: gf      !! g value at tf

    !local variables:
    real(wp) :: t,dt,t2,ga,gb,dt_root,dum
    real(wp),dimension(me%n) :: x !! state vector
    real(wp),dimension(me%n) :: g_xf !! state vector from the root finder
    logical :: first !! it is the first step
    logical :: last  !! it is the last step
    integer :: iflag !! return flag from `solver`

    if (.not. associated(me%f)) then
        call me%raise_exception(RKLIB_ERROR_F_NOT_ASSOCIATED)
        return
    end if
    if (.not. associated(me%g)) then
        call me%raise_exception(RKLIB_ERROR_G_NOT_ASSOCIATED)
        return
    end if
    if (abs(h)<=zero) then
        call me%raise_exception(RKLIB_ERROR_INVALID_H)
        return
    end if

    call me%begin_integration()

    call me%export_point(t0,x0,.true.) !first point

    if (abs(t0-tmax)<=zero) then
        xf = x0
        tf = t0
        call me%g(t0,x0,gf)
    else

        first = .true.
        t = t0
        x = x0
        call me%g(t0,x0,ga)     !evaluate event function
        dt = sign(h,tmax-t0)    !time step (correct sign)

        do

            t2 = t + dt
            last = ((dt>=zero .and. t2>=tmax) .or. &  !adjust last time step
                    (dt<zero .and. t2<=tmax))         !
            if (last) then
                dt = tmax-t
                t2 = tmax
            end if
            call me%step(t,x,dt,xf)
            if (me%stopped) return
            me%num_steps = me%num_steps + 1
            if (me%num_steps > me%max_number_of_steps) then
                call me%raise_exception(RKLIB_ERROR_TOO_MANY_STEPS)
                return
            end if
            call me%g(t2,xf,gb)     !evaluate event function

            if (first .and. abs(ga)<=tol) then

                !we ignore a root at t0 after the first step
                if (abs(gb)<=tol) then !check this one since it could have landed on a root
                    gf = gb
                    tf = t2
                    exit
                else
                    if (last) then  !exiting without having found a root
                        tf = t2
                        gf = gb
                        exit
                    end if
                    call me%export_point(t2,xf)   !intermediate point
                    x = xf
                    t = t2
                    ga = gb
                end if

            elseif (ga*gb<=zero) then !there is a root somewhere on [t,t+dt]

                !find the root:
                call root_scalar(me%solver,solver_func,zero,dt,dt_root,dum,iflag,&
                                 fax=ga, fbx=gb, rtol=tol, atol=tol)
                                ! ftol,maxiter,bisect_on_failure) ! other options
                if (me%stopped) return
                t2 = t + dt_root
                gf = solver_func(dt_root)
                if (me%stopped) return
                tf = t2
                xf = g_xf !computed in the solver function
                exit

            else  !no root yet, continue

                if (last) then  !exiting without having found a root
                    tf = t2
                    gf = gb
                    exit
                end if
                call me%export_point(t2,xf)   !intermediate point
                x = xf
                t = t2
                ga = gb

            end if

            if (first) first = .false.

        end do

    end if

    call me%export_point(t2,xf,.true.)   !last point

    contains

        function solver_func(delt) result(g)

        !! root solver function. The input is the `dt` offset from time `t`.

        implicit none

        real(wp),intent(in) :: delt  !! from [0 to `dt`]
        real(wp) :: g

        !take a step from t to t+delt and evaluate g function:
        call me%step(t,x,delt,g_xf)
        if (me%stopped) return
        call me%g(t+delt,g_xf,g)

        end function solver_func

    end subroutine integrate_to_event_fixed_step