Nonlinear ODE system solver (Runge-Kutta)

Question

I`m trying to solve the following nonlinear-system(please see the pic attached) but I´m truly lost.

EQUATIONS

I need to create one method called solve_prec, which solves the system. I dont know how to use scipy.integrate.solve_ivp in this case.

Note: T (x=0) = 1, and the domain of integration is (0, -0.25). Moreover energy method refers to e, pressure method refers to p. Finally $\rho$ is den, e0 = energy(1,1) and p0 = pressure(1,1) Thanks in advance!!

class eos_gamma_c():

 def __init__(self, gamma=5/3.0, M0 = 10, P0 = 0.0001, k = 0.0001):
  self.gamma = gamma
  self.M0 = M0
  self.P0 = P0
  self.k  = k
  self.cv = 1./self.gamma/(self.gamma-1.)
  return

def energy(self, tem, den):
 return self.cv*tem 

def pressure(self, tem, den):
 return den*tem/self.gamma

def solve_prec(self, tem, den):
 ene0 = self.energy(1,1)
 pre0 = self.pressure(1,1)