compute_stepsize Subroutine

    subroutine compute_stepsize(me,n,h,tol,err,p,hnew,accept)

    implicit none

    class(stepsize_class),intent(inout) :: me  !! (`inout` since it has some state, see `pi_beta`
                                               !! and `no_increase_after_reject`)
    integer,intent(in)               :: n      !! number of variables
    real(wp),intent(in)              :: h      !! current step size (<>0)
    real(wp),dimension(n),intent(in) :: tol    !! abs error tolerance (>0)
    real(wp),dimension(n),intent(in) :: err    !! truncation error estimate (>0)
    integer,intent(in)               :: p      !! order of the method
    real(wp),intent(out)             :: hnew   !! new step size (<>0)
    logical,intent(out)              :: accept !! if the step is accepted

    real(wp) :: e        !! exponent
    real(wp) :: hfactor  !! step size factor (>0)
    real(wp) :: max_err  !! max error for all the elements

    real(wp),parameter :: small = 10.0_wp * epsilon(1.0_wp) !! small value of the scaled error

    if (me%fixed_step_mode) then
        ! do not adjust the step size
        accept = .true.
        hnew = h
    else

        ! the scaled error (the step is acceptable if <= 1):
        if (me%relative_err) then
            max_err = me%norm(err/(tol*abs(h)))  ! error per unit step
        else
            max_err = me%norm(err/tol)           ! error per step
        end if

        if (max_err <= small) then ! the error is extremely small (relative to the tolerance)

            hfactor = me%hfactor_accept
            accept = .true.

        else

            e = 1.0_wp / real(max(1, p+me%p_exponent_offset),wp)
            e = e - me%pi_beta_mult*me%pi_beta  ! (if using the PI control)
            hfactor = abs( me%safety_factor * (1.0_wp/max_err)**e )

            ! if the step is to be accepted:
            select case (me%accept_mode)
            case(1) !algorithm 17.12
                accept = hfactor >= 1.0_wp
            case(2) !algorithm 17.13
                accept = max_err <= 1.0_wp
            end select

            !...notes:
            ! see: L. Shampine "Some Practical Runge-Kutta Formulas",
            !      Mathematics of Computation, 46(173), Jan 1986.
            ! different conditions for satisfying error conditions:
            !  ||err|| <= tol   -- Error per step (EPS)
            !  ||err|| <= h*tol -- Error per unit step (EPUS)

            !compute the actual hfactor based on the limits:
            if (accept) then
                ! PI control: include the error of the previous accepted step
                if (me%pi_beta > zero) hfactor = hfactor * me%err_old**me%pi_beta
                hfactor = min(me%hfactor_accept, hfactor)
            else
                hfactor = max(me%hfactor_reject, hfactor)
            end if

        end if

        if (accept) then
            if (me%no_increase_after_reject .and. me%last_rejected) hfactor = min(hfactor, 1.0_wp)
            me%err_old = max(max_err, 1.0e-4_wp)
        end if
        me%last_rejected = .not. accept

        ! compute the new step size (enforce min/max bounds & add sign):
        hnew = sign(max(me%hmin,min(me%hmax,abs(h)*hfactor)),h)

    end if

    end subroutine compute_stepsize