Hello, I am trying to now use my Runge Kutta code to numerically solve the equations of state for an ordinary star and then transcribe that to a neutron star. My code is below: Code (Text): import numpy as np import matplotlib.pyplot as plt from math import pi #Constants G = 6.674*1e-11 N = 5000 Gamma = 1.5 K = 5000 dur = 5 dr = dur/float(N-1) state = np.zeros([N,2]) #initial conditions state[0,0] = 0 state[0,1] = 5 def rk4(y,r, dr, den, deriv): k1 = dr * deriv(y, den, r) k2 = dr * deriv(y + 0.5*k1,den, r + 0.5*dr) k3 = dr * deriv(y + 0.5*k2,den, r + 0.5*dr) k4 = dr * deriv(y + k3, den, r + dr) y_next = y + (k1 + 2*(k2+k3)+k4)/6 return y_next def density(state): return ((1/K)*state)**(-Gamma) # row = (P/K)^(-Gamma) def EoS(state,den, r): y0 = 4*pi*(r**2)*den #dM/dr = 4*pi*(r^2)*row y1 = (-G*state*den)/(r**2) #dP/dr = -G*M*row / r^2 return np.array([y0,y1]) r = 1e-23 #starting value to avoid singularity for i in range(N-1): den = density(state[i]) state[i+1]= rk4(state[i],r,dr,den, EoS) r = r + dr plt. plot(state[:,1], state[:,0]) plt. show() What I am trying to do is solve the two equations of state while using the polytropic relation between density and pressure.Now, I am slightly confused on how I can accomplish this as I have three unknowns: mass, pressure, and density which requires three separate equations are needed to develop a unique solution. I, however, am only receiving a straight line that diverges to infinity on the second iteration. Does anyone have any ideas why this is occurring?