lowable step size.
f0 : ndarray, shape (n,)
    Initial value of the derivative, i.e., ``fun(t0, y0)``.
direction : float
    Integration direction.
order : float
    Error estimator order. It means that the error controlled by the
    algorithm is proportional to ``step_size ** (order + 1)`.
rtol : float
    Desired relative tolerance.
atol : float
    Desired absolute tolerance.

Returns
-------
h_abs : float
    Absolute value of the suggested initial step.

References
----------
.. [1] E. Hairer, S. P. Norsett G. Wanner, "Solving Ordinary Differential
       Equations I: Nonstiff Problems", Sec. II.4.
r