Note that the code I am using is Python. Currently, I am attempting to numerically solve the equations of state that give rise to the Chandrasekar limit for a white dwarf star. My code works for a simple harmonic oscillator given the correct equations of motion, but does not compute the correct curve for my two coupled ODEs that I need for my project. My code is below. Note that the equations of state are defined in the comments. Code (Python): import numpy as np import matplotlib.pyplot as plt from math import sqrt #Variable Definitions N = 500 #number of entries dur = 25000 #radius taken up to in kilometers dr = dur / float(N-1) #step size #creating the array y = np.zeros([N,2]) y[0,0] = 20 #initial central momentum y[0,1] = 0 #initial mass #runge kutta algorithm def rk4(y, r, dr, deriv): k1 = dr * deriv(y,r) k2 = dr * deriv(y + 0.5*k1,r + 0.5*dr) k3 = dr * deriv(y + 0.5*k2, r + 0.5*dr) k4 = dr * deriv(y + k3, r + dr) y_next = y + (k1 + 2*(k2+k3)+k4)/6 return y_next #Equations of state (Note that x is my fermi momentum x and x is my mass. r is my radius def Eqns(x,r): y0 = (-5*x*sqrt(1+x**2))/(3*x*r**2) #dx/dr = -5/3(M/r^2)*sqrt(1+x^2)/x y1 = 3*(r**2)*(x**3) #dM/dr = 3*r^2*x^3 return np.array([y0,y1]) #forming data points for i in range(N-1): y[i+1] = rk4(y[i],10**(-23), dr, Eqns) #started at some arbitrary smalll value to avoid #singularity radius = np.linspace(10**(-23),25000, N) #plotting plt.plot(y[:,1],radius) plt.show() For some reason when I run this all I receive is a straight linear line which is nothing like the curve I expected that is attached to this post. Also, the graph indicates my mass is negative which is ridiculous. I have a feeling the error is a semantic error somewhere in my code, but would I not expect a curve instead of a straight line for ODEs that are obviously not linear. Thanks again. EDIT: I forgot to mention that what the curve should resemble is the curve in the relativistic fermi gas line.