Hi! Thanks in advance for your help!

I'm currently working in my dissertation on Neutron Stars in f(R) gravity. Therefore, as a first approach, I wanted to solve TOV equations in general relativity using the usual formulation with {p', m'} to get M-R plots. I'm well aware of the singularity at r = 0 for the equation for p, hence, I'm trying to implement it from a r = rmin value rather than r = 0. Nevertheless, I keep getting "Step size is effectively zero; singularity or stiff system suspected." and curves different than the ones expected.

I've already tried using Piecewise in order to modify the p' equation or adding a small amount to r, as r = r + 0.001, for example, but the result isn't the expected.   I've seen this being solved in Python and C using this same approach, with a rmin value for r with good results and plots so I expected this to be easily implementable in Mathematica as well, but I haven't been able to. Any help would be highly appreciated!

Expected result or similar. This wasn't get using the usual formalism but rather Lindblom formalism where an enthalpy variable is introduced (Physics stack exchange forum).