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.
