  In this paper, we study the thermalization of gluons and N_f flavors of
massless quarks and antiquarks in a spatially homogeneous system. First, two
coupled transport equations for gluons and quarks (and antiquarks) are derived
within the diffusion approximation of the Boltzmann equation, with only 2<-> 2
processes included in the collision term. Then, these transport equations are
solved numerically in order to study the thermalization of the quark-gluon
plasma. At initial time, we assume that no quarks or antiquarks are present and
we choose the gluon distribution in the form f = f_0 theta (1-p/Q_s) with Q_s
the saturation momentum and f_0 a constant. The subsequent evolution of systems
may, or may not, lead to the formation of a (transient) Bose condensate,
depending on the value of f_0. In fact, we observe, depending on the value of
f_0, three different patterns: (a) thermalization without gluon Bose-Einstein
condensate (BEC) for f_0 < f_{0t}, (b) thermalization with transient BEC for
f_{0t} < f_0 < f_{0c}, and (c) thermalization with BEC for f_{0c} < f_0. The
values of f_{0t} and f_{0c} depend on N_f. When f_0> 1 > f_{0c}, the onset of
BEC occurs at a finite time t_c ~ 1/((alpha_s f_0)^2 Q_s). We also find that
quark production slows down the thermalization process: the equilibration time
for N_f = 3 is typically about 5 to 6 times longer than that for N_f = 0 at the
same Q_s.
