  We study the Wouthuysen-Field coupling at early universe with numerical
solutions of the integrodifferential equation describing the kinetics of
photons undergoing resonant scattering. The numerical solver is developed based
on the weighted essentially non-oscillatory (WENO) scheme for the
Boltzmann-like integrodifferential equation. We focus on the time evolution of
the Wouthuysen-Field (W-F) coupling in relation to the 21 cm emission and
absorption at the epoch of reionization. We show that a local Boltzmann
distribution will be formed if photons with frequency \sim \nu_0 have undergone
a ten thousand or more times of scattering, which corresponds to the order of
10^3 yrs for neutral hydrogen density of the concordance \Lambda CDM model. The
time evolution of the shape and width of the local Boltzmann distribution
actually doesn't dependent on the details of atomic recoil, photon sources, or
initial conditions very much. However, the intensity of photon flux at the
local Boltzmann distribution is substantially time-dependent. The time scale of
approaching the saturated intensity can be as long as 10^5-10^6 yrs for typical
parameters of the \Lambda CDM model. The intensity of the local Boltzmann
distribution at time less than 10^5 yrs is significantly lower than that of the
saturation state. Therefore, it may not be always reasonable to assume that the
deviation of the spin temperature of 21 cm energy states from cosmic background
temperature is mainly due to the W-F coupling if first stars or their
emission/absorption regions evolved with a time scale equal to or less than
Myrs.
