How do I differentiate this complicated function? (The Poisson equation for the isothermal sphere)?

Question

Ignore the first part of the question. Look at part ii. I understand how to solve basic differential equations in python, but when you have a nested differential I get quite confused. Please help

Answer 1

I will use a different first-order system using the nested structure of the given equation, set z=x^2*d(ln y)/dx, which implies z(0)=0 for the assumed-as-non-singular solution y, then

y'=z*y/x^2
z'=-6*x^2*y

At x close to zero we have by the initial condition y close to 1 so that the second equation integrates in the lowest degree terms to z=-2*x^3. Inserted into the first equation this gives indeed y=1-x^2. One could now continue to compute higher-degree terms...

Now use the standard tools of python-scipy as demonstrated in the documentation or any other source for, among others, the popular Lorenz model as example of a multi-dimensional coupled system, leading to a code of

x = np.logspace(-1,3,500)
x0 = x[0]
y0 = 1-x0**2
z0 = -2*x0**3

def odesys(x,u):
    y,z = u
    return z*y/x**2, -6*x**2*y

u = odeint(odesys, [y0,z0], x, tfirst=True, atol=1e-9, rtol=1e-11)
y,z = u.T

Putting the computed solution into a loglog plot along with a plot of x^-2 gives the graphs

which shows the claimed trend.