(*
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