Solving Ode With Scipy Gives Invalid Value Encountered In Double_scalars Error

I'm trying to solve the Lane-Emden equation for arbitrary values of n (polytropic index). In order to use SciPy, I have represented the second-order ODE as a set of two coupled fir

Solution 1:

The following works for me (which looks similar to the Wikipedia entry). The derivative was written incorrectly (negative on the wrong element) and input to the derivative is changed simply to n.

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint

defpoly(y,xi,n):  
    theta, phi = y  
    dydxi = [-phi/(xi**2), (xi**2)*theta**n]  
    return dydxi

fig, ax = plt.subplots()

y0 = [1.,0.]  
xi = np.linspace(10e-4, 16., 201)

for n inrange(6):
    sol = odeint(poly, y0, xi, args=(n,))
    ax.plot(xi, sol[:, 0])

ax.axhline(y = 0.0, linestyle = '--', color = 'k')
ax.set_xlim([0, 10])
ax.set_ylim([-0.5, 1.0])
ax.set_xlabel(r"$\xi$")
ax.set_ylabel(r"$\theta$")
plt.show()

---

The original question used fractional values of the polytropic index, so something like follow could be used to investigate those values ...

for n in range(12):
    sol = odeint(poly, y0, xi, args=(n/2.,))

    if np.isnan(sol[:,1]).any():
        stopper = np.where(np.isnan(sol[:,1]))[0][0]
        ax.plot(xi[:stopper], sol[:stopper, 0])
    else:
        ax.plot(xi, sol[:, 0])