The temporal relaxation of the electron component in weakly ionized, collision-dominated plasmas is studied on the basis of the time-dependent Boltzmann equation. The relaxation process is treated using an expansion of the electron velocity distribution function in Legendre polynomials. In the past this expansion has been truncated after the first two terms. The equation system for the first two expansion coefficients which is obtained in this two-term approximation from the time-dependent kinetic equation has been solved by the so-called conventional approach and, very recently, by a strict time-dependent approach. In the conventional approach the distribution anisotropy is treated in a quasi-steady-state way. In generalization of the conventional and the strict time-dependent two-term approximation an efficient numerical method for solving the time-dependent electron kinetic equation in a multi-term approximation has been developed and its main aspects are presented in this paper. As an application, the temporal relaxation of an energetic electron group in a weakly ionized, collision-dominated neon plasma, acted upon by a DC electric field, has been investigated. Main findings of this relaxation study are that, depending on the reduced electric field strength, (i) the conventional two-term approximation can completely break down at early and intermediate stages of the relaxation process, (ii) the strict time-dependent two-term approximation can also be insufficient at intermediate relaxation stages and (iii) the converged solution of the Boltzmann equation is obtained throughout the relaxation in an approximation with six to eight terms from the expansion of the electron velocity distribution.