Scipy solution of coupled ODEs flatlines after a critical value

I am trying to solve a system of two coupled ODEs using Scipy's solve_ivp function. Namely the hydro-static equilibrium and the mass for a white dwarf with full units.

My initial conditions are that the enclosed mass should be 0 at the core, and that the pressure should be 1e24 at the center. The code works for central pressures up to 1e16 but past that critical point, the solution for the pressure flatlines. The initial value does not change.

When I compile the code, there are no errors nor any warnings produced.

import numpy as np
import math as m
import matplotlib.pyplot as plt

from scipy.integrate import solve_ivp

def dSdx(r, S):
    p, m = S
    mu = 2 # This is an approx. 
    gamma = 4/3
    K = 1.2435e15/(mu**gamma)
    G = 6.67430e-8 # cgs [cm^3 g^-1 s^-2]
    if r==0:
        return [0,0]
    else:
        return [-m*G* K**(gamma)/(r**2 * p**gamma), 4* K**(gamma)* np.pi *r**2 /p**gamma]

S_0 = [1e24,0]
dt=10000
sol = solve_ivp(fun=dSdx, t_span=(0,1e9), max_step=dt, y0=S_0, method='RK45', 
                atol=0.01, rtol=0.001)

p_sol = sol.y[0] 
m_sol = sol.y[1]/(1.989*1e33) # Solar Masses
x=sol.t*1e-5 # Km

# Plotting
plt.figure(1)
fig = plt.plot(x,p_sol)
plt.legend(['Pressure'])
plt.title('Pressure of a Newtonian White Dwarf')
plt.xlabel("R $[Km]$")
plt.ylabel("Pressure $[erg/cm^3] $")

plt.figure(2)
fig = plt.plot(x,m_sol, color='darkorange')
plt.legend(['Mass'])
plt.title('Mass of a Newtonian White Dwarf')
plt.xlabel("R $[Km]$")
plt.ylabel("Mass $[M_{\oplus}]$")