  The study of strongly out-of-equilibrium states and their time evolution
towards thermalization is critical to the understanding of an ever widening
range of physical processes. We present a numerical method that for the first
time allows for the numerical solution of the most difficult part of the
time-dependent Boltzmann equation: the full scattering term. Any number of
bands (and quasiparticles) with arbitrary dispersion, any number of high order
scattering channels (we show here four legs scatterings: electron-electron
scattering) can be treated far from equilibrium. No assumptions are done on the
population and all the Pauli-blocking factors are included in the phase-space
term of the scattering. The method can be straightforwardly interfaced to a
deterministic solver for the transport. Finally and most critically, the method
conserves to machine precision the particle number, momentum and energy for any
resolution, making the computation of the time evolution till complete
thermalization possible. We apply this approach to two examples, a metal and a
semiconductor, undergoing thermalization from a strongly out-of-equilibrium
laser excitation. These two cases, which are in literature treated hitherto
under a number of approximations, can be addressed free from those
approximations and straightforwardly within the same numerical method.
