Solving the TOV equations does indeed require simultaneously solving two ODEs:
\begin{align}
\frac{\mathrm{d}P}{\mathrm{d}r}&=TOV(P,\,\rho,\,m) \\
\frac{\mathrm{d}m}{\mathrm{d}r}&=4\pi\rho r^2
\end{align}
However, this cannot be solved as-is because you have two equations and three unknowns (mass, pressure & density). Hence, you also need to make an assumption on the functional form of the equation of state (EOS), $P(\rho)$ which, assuming you have a nice, analytic function, can be invertible to find the density from the pressure (i.e., $\rho(P)$). A common function, for instance, is the polytrope ($P(\rho)=\kappa\rho^\gamma\leftrightarrow\rho=(P/\kappa)^{1/\gamma}$).
If your EOS is not nice, you may require finding the density from the pressure via root-finding (i.e., finding $\hat{\rho}$ s.t. $P_\text{input}-P(\hat{\rho})\simeq0$).
You then solve the system much the same way that FlatterMann describes. A somewhat simplified example in Python is below.
def eos(rho: float) -> float:
""" equation of state """
return k * pow(rho, gamma)
def inv_eos(p: float) -> float:
""" inverse equation of state """
return pow(p / k, 1 / gamma)
def tov_integrate(y: ArrayLike, r: float) -> ArrayLike:
""" integration of the TOV equation """
p, m = y[0], y[1]
rho = inv_eos(p)
dp_dr = TOV(p, m, rho)
dm_dr = 4 * pi * r*2 rho
return [dp_dr, dm_dr]
def integrate(rho_central: float) -> tuple
""" integrate the system """
# set up r as the range [?, R] with N bins
m = zeros(len(r))
p = zeros(len(r))
m[0] = 4 * pi * r[0]*2 rho_central
p[0] = eos(rho_central)
i, y = 0, [p[0], m[0]]
while p[i] > 0 and i < N:
y = rk4(tov_integrate, y, r[i], r[i+1] - r[i])
p[i+1] = y[0]
m[i+1] = y[1]
i += 1
return r, p, m
You may also want to save and return the densities for plotting purposes, but that could be determined using the EOS anyway.
