Using a solitonic model of nuclear matter, the Bogomol'nyi-Prasad-Sommerfield (BPS) Skyrme model, we compare neutron stars obtained in the full field theory, where gravitational backreaction is completely taken into account, with calculations in a mean-field approximation using the Tolman-Oppenheimer-Volkoff approach. In the latter case, a mean-field-theory equation of state is derived from the original BPS field theory. We show that in the full field theory, where the energy density is nonconstant even at equilibrium, there is no universal and coordinate-independent equation of state of nuclear matter, in contrast to the mean-field approximation. We also study how neutron star properties are modified by going beyond mean-field theory and find that the differences between mean-field theory and exact results can be considerable. Further, we compare both exact and mean-field results with some theoretical and phenomenological constraints on neutron star properties, demonstrating thus the relevance of our model even in its most simple version.