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)