I'm trying to write a Python code to implement Hartle-Thorne approximation formalism for rotating neutron stars. In this formalism, I am required to solve the following second order differential equation: 1r4ddr(r4dω(r)ˉdr)+4rdj(r)drω(r)ˉ=0 \frac{1}{r^4}\frac{d}{dr}(r^4 \frac{d \bar{\omega(r)}}{dr}) + \frac{4}{r}\frac{dj(r)}{dr}\bar{\omega(r)} = 0 with the boundary conditions ωˉ(r=0)=ωcˉ\bar{\omega}(r=0) = \bar{\omega_c} and dωˉdr∣r=0=0\frac{d\bar{\omega}}{dr}|_{r=0} = 0 where ωcˉ\bar{\omega_c} is an arbitrary constant. jj is defined as j=e−ν(r)1−2m(r)rj = e^{-\nu(r)}\sqrt{1-\frac{2m(r)}{r}} inside the star and j(r)=1j(r)=1 for r≥Rr \geq R where RR is the radius of the star . ν\nu is the metric function defined as dνdr=−2(ϵ+p)−1dpdr\frac{d\nu}{dr} = -2(\epsilon + p)^{-1}\frac{dp}{dr} I obtained the metric by simultaneously solving it's differential equation with the TOV equations using a RK4 solver. Now I am trying to solve that second order differential equation using odeint but the result does not seem to be correct. The issues are as follows: The results I obtained does not match with the results published in the related paper. The results change if I change the boundary condition on ωˉ\bar{\omega} even though this boundary condition is supposed to be arbitrary. I need some help to identify the correct approach to solve this equation: can I solve it separately now that I already have obtained all the necessary quantities from the TOV-RK4 run (that's what I am doing), or do I need to solve it with the TOV equations themselves? For those who might need it, here is the paper that describes the formalism and all the necessary equations.