def rk4_step(state, time, tau, state_derivative_function, **kwargs):
    """
    Routine for Rungge-Kutta steps.
    """
    fun = state_derivative_function
    k1 = tau * fun(state, time, **kwargs)
    state_ = state + 0.5 * k1
    k2 = tau * fun(state_, time + tau * 0.5, **kwargs)
    state_ = state + 0.5 * k2
    k3 = tau * fun(state_, time + tau * 0.5, **kwargs)
    state_ = state + k3
    k4 = tau * fun(state_, time + tau, **kwargs)
    ans = state + k1 / 6 + k2 / 3 + k3 / 3 + k4 / 6
    return ans
