What am I doing wrong with my multidimensional state array?

[Bishop556]

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:

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[1])**(-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[1]*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 + drplt. 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?

 

[Bishop556]

Also, note that state[0] is mass, state[1] is pressure, and den is density.

 

[Mark44]

Also, note that state[0] is mass, state[1] is pressure, and den is density.

I'm not sure about this, but it seems to me that you defined your state array two be two-dimensional (your code: state = np.zeros([N,2]) ), but you are using it later as if it were a one-dimensional array, as in what you wrote above.