The BPS Skyrme model has been demonstrated already to provide a physically intriguing and quantitatively reliable description of nuclear matter. Indeed, the model has both the symmetries and the energy–momentum tensor of a perfect fluid, and thus represents a field theoretic realization of the “liquid droplet” model of nuclear matter. In addition, the classical soliton solutions together with some obvious corrections (spin–isospin quantization, Coulomb energy, proton–neutron mass difference) provide an accurate modeling of nuclear binding energies for heavier nuclei. These results lead to the rather natural proposal to try to describe also neutron stars by the BPS Skyrme model coupled to gravity. We find that the resulting self-gravitating BPS Skyrmions provide excellent results as well as some new perspectives for the description of bulk properties of neutron stars when the parameter values of the model are extracted from nuclear physics. Specifically, the maximum possible mass of a neutron star before black-hole formation sets in is a few solar masses, the precise value of which depends on the precise values of the model parameters, and the resulting neutron star radius is of the order of 10 km.