def _udot(u, t, coeff_rec_func, redshift_func, ion_func, bubble=True):
    """du/dt where u = x - f_* f_esc,gamma N_gamma F
    
    Parameters
    ----------

    u: integral of du/dt as defined below

    t: cosmic age in s

    redshift_func: function returning redshift given t
    
    ion_func: function returing ionized fraction neglecting recombinations

    coeff_rec_func: function returning clumping_factor alpha_B n_H_0 (1+z)^3

    bubble: If True, assume ionized gas is in fully-ionized bubbles
            and other gas is fully neutral. If False, asssume gas is
            uniformly fractionally ionized.

    Notes
    -----

    This is implemented as a reformulation of the normal ODE
    describing ionization and recombination (see, e.g. Bagla, Kulkarni
    & Padmanabhan (2009MNRAS.397..971B).

    The original ODE is:

    dx/dt = -alpha_B C n_H x + f_* f_esc,gamma N_gamma dF/dt

    If we let u = x - w, where w = f_* f_esc,gamma N_gamma F(t) then

    du/dt = dx/dt - dw/dt

    which gives

    du/dt = -alpha_B C n_H x = -alpha_B C n_H (u + w)

    We have an analytical expression for w, so we can numerically
    integrate the ODE to give us u(t) or x(t) = u(t) + w(t).

    """
    z = redshift_func(t)
    w = ion_func(z)
    crf = coeff_rec_func(z)
    x = u + w
    x = x * (x <= 1.0) + 1.0 * (x > 1.0)
    if bubble:
        udot = -1.0 * crf * x
    else:
        udot = -1.0 * crf * x ** 2
    if False:
        print('z=%.3f; t=%.1g; c=%.2g; udot=%.2g; w,x,u = %.2g, %.2g, %.2g' %
            (z, t, crf, udot, w, x, u))
    return udot
