# Numerical solution to the Tolman-Oppenheimer-Volkoff equations for any equation of state (numerical or analytical)

I've been working on a code to solve the Tolman-Oppenheimer-Volkoff (TOV) equations for a while and recently I've got it right but only for one specific equation of state, the bag model, which is not what I need since it should run for any hadronic equation of state. The code has given only wrong results for any EoS other than the bag model and I've tried everything I could think of but nothing has made it so that the final diagrams for other EoS are correct. I'll put the code and some of the wrong results for some of the EoS I've been studying below.

import math
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.interpolate import UnivariateSpline
from scipy.interpolate import interp1d

pi = math.pi
km2_gevfm3 = (1/(1.3234*10**(-6)))*0.001
gevfm3_km = 1/(km2_gevfm3)
B = 7.424*10**(-5)
cm = 1/1.4760

eos = np.genfromtxt("C:\\Users\\felip\\Desktop\\TOV\\Equações de estado\\PtableTM1.dat")
eden, p = (gevfm3_km*eos[:,0])/1000, (gevfm3_km*eos[:,1])/1000
press = interp1d(eden, p, kind = 'cubic', bounds_error = False, fill_value = 'extrapolate')
# press = UnivariateSpline(eden, p, k = 5, s = 5)
dpde = press._spline.derivative(nu = 1)
# dpde = press.derivative(n = 1)

def TOV(y, r):
eden, M = float(y[0]), float(y[1])
dedr = -(1*(dpde(eden))**(-1))*((M*eden)/(r**2))*(1+press(eden)/eden)*(1+(4*pi*press(eden)*r**(3))/(M))*((1-((2*M)/r))**(-1))
dMdr = 4*pi*eden*(r**2)
return [dedr, dMdr]

r0 = 10**(-10)
rm = 5000  #km
N = 1000
rn = np.linspace(r0, rm, N)

emin, emax =   7.7*10e-15, 0.0003
h_e = 1000
e_range = np.linspace(emin, emax, h_e)
R,M,E =[],[],[]

for e0 in e_range:
M0 = (4*pi*e0*r0**(3))/3
contorno = [e0, M0]
soltov = odeint(TOV, contorno, rn, atol = 10e-10, rtol = 10e-14, mxstep=5000)
rmax = r0
mmax = M0
esol = soltov[:,0]
psol = press(esol)
k=0
for i in psol[:-1]:
h=rn[k+1]-rn[k] + r0
rmax = rmax + h
mmax = soltov[k,1]
k=k+1
if i <= 0:
break
R.append(rmax)
M.append(mmax*cm)
E.append(e0*km2_gevfm3)

immax = np.argmax(M)
Mmax = M[immax]
Rmax = R[immax]

f1 = plt.figure()
ax1.title.set_text('M x E')
ax1.set_xlabel('$$(\epsilon)_c [GeV/fm^3]$$')
ax1.set_ylabel('$$M [M_\odot]$$')
ax1.plot(E,M,'-',color='blue')
ax1.legend()

f2 = plt.figure()
ax2.set_xlabel('$$R [Km]$$')
ax2.set_ylabel('$$M [M_\odot]$$')