Please see the attached plot of this type for your reference. How to plot M-R curves by solving a coupled ordinary differential equations using "NDSolve" package?
Mass - Radius relation for compact object provides the variation of mass w.r.t. to its radius which finally gives the unique upper limit or maximum mass . For example , compact object white dwarf has maximum mass limit is 1.4 Mo (Mo being the solar mass ). For Neutron Star we do not have any upper limit or maximum mass. The reason is that we do not know what is actually there in the interior as well as at the center of the neutron star (?). Interior of neutron star contains superconductor, superfluid form of protons, neutrons etc that are almost certain but exactly not yet certain. At the center it is suggested that core consists of quark matter . In other words , neutron star with quark matter core type behaviour has recently been observed by Anna et al published in nature, 2020. It is also suggested that under extreme pressure and density the deconfined phase i.e. u (up), d( down), s(strange) quarks are present with asymptotic freedom phase i.e. u, d, s quarks are roaming freely and strong forces are vanished.
Different models have been proposed based on various constituents and their interactions. Naturally, we shall get various maximum mass limit by plotting M-R for different assumptions regarding the constituents in the interior of neutron star. Until and unless we will certain about the constituents at the center of neutron star. We have to think an alternative way to find out the maximum mass of neutron star by avoiding the interactions. I have proposed an alternative way and estimated maximum mass of neutron star is = 3.199 Mo i.e. maximum mass of neutron star can not more than 3.2 Mo.
A straightforward answer is: use NDSolve to solve the TOV equations setting the boundary conditions at a very small radial coordinate. Use the option "WhenEvent" to stop the integration when the pressure reaches zero.