The most probable velocity distribution function of each component,fa, of a nonequilibrium multispecies spherically symmetric system of particles (stellar plasma atmospheres and winds, stellar systems, pellet-fusion systems) is analytically derived forthe case in which each component is described by the first six moments of fa. This is achieved by the aid of a variational approach based on the requirement that the BoltzmannH function for the system be a minimum, subject to the constraints provided by the sets of six macroscopic parameters describing the nonequilibrium state. The use of the so-obtained velocity distribution functions for the closure of the moment equations as well as for the calculation of their collisional terms (via the Fokker-Planck equation) is discussed. The limitations on the maximum deviations from the equilibrium state which are consistent with the assumptions used are also indicated.