(* Solution to Problem 2 of the Bigger Problem Set - Lane-Emden equation. At the Mathematica prompt, type << laneemden.math *) (* The polytropic index to do numerical calculations for *) nn = 5; (* maximum x to integrate to *) xp = 100; (* Series solution *) s[n_,k_,x_] := Sum[a[n,i] x^i , {i,0,k}]; (* Coefficients of series solution *) a[n_, k_ /; k > 0] := a[n,k] = \ - Coefficient[ Series[s[n,k-2,x]^n, {x,0,k-2}], x, k-2] / (k (k+1)); (* Series solution in the particular case a_0 = 1 *) ths[n_,k_,x_] := s[n,k,x] /. a[n,0] -> 1 (* The Lane-Emden equation, with boundary condition from the series solution evaluated at x = e *) deq[n_,e_] := {th''[x] + 2/x th'[x] + th[x]^n == 0, th[e] == ths[n,12,e], th'[e] == D[ths[n,12,x],x] /. x -> e}; (* Solve the diffeq for the case n = nn starting at small x *) xm = .01; soln = NDSolve[deq[nn,xm], th, {x,xm,xp}]; (* Plot the solution th, and the dimensionless density th^n *) Plot[{th[x] /. soln, th[x]^nn /. soln}, {x,xm,xp}]; (* Value of x where th = 0 *) xedge = x /. FindRoot[ Re[th[x] /. soln][[1]], {x,{xm,xp}}] (* Derivative of th at the edge *) thpedge = Re[ D[th[x] /. soln, x] /. x -> xedge ] (* Ratio rhobar/rhoc of mean to central density *) rhobartorhoc = - 3 thpedge/xedge