This paper is devoted to multiplicity results of solutions to nonlocal
elliptic equations modeling gravitating systems. By considering the case of
Fermi-Dirac statistics as a singular perturbation of Maxwell-Boltzmann one, we
are able to produce multiplicity results. Our method is based on cumulated mass
densities and a logarithmic change of coordinates that allows us to describe
the set of all solutions by a non-autonomous perturbation of an autonomous
dynamical system. This has interesting consequences in terms of bifurcation
diagrams, which are illustrated by a some numerical computations. More
specifically, we study a model based on the Fermi function as well as a
simplified one for which estimates are easier to establish. The main difficulty
comes from the fact that the mass enters in the equation as a parameter which
makes the whole problem non-local
