(* The solution is obtained by integrating Lane-Emden, The intigration begins form ξ=ϵ because the prefactor vanishes at ξ=0 *)

ϕSol[n_]:=
Module[{ξ,ϕ,ϵ=$MachineEpsilon},
NDSolve[{ξ ((ϕ[ξ])^n)+2 ϕ'[ξ]+ξ ϕ''[ξ]==0,ϕ[ϵ]==1,ϕ'[ϵ]==0},ϕ,{ξ,ϵ,2+n*3}][[1,1,2]]
]

(* This is the solution for n=1.5 and for n=3 *)

Plot[{Evaluate[ϕSol[1.5][ξ]],Evaluate[ϕSol[3][ξ]]},{ξ,0,8},Frame->True,PlotRange->{Automatic,{0,1}},FrameLabel->{"ξ","ϕ"},LabelStyle->{Directive[FontSize->18],Directive[FontSize->18]},FrameStyle->{{Thickness[0.001],Thickness[0.001]},{Thickness[0.001],Thickness[0.001]}},Frame->True,AspectRatio->1/1.618]

(* This is the first zero for a polytrope with index n *)

ξ1[n_]:=Abs[ξ/.FindRoot[ 0==Evaluate[ϕSol[n][ξ]],{ξ,3}]];

(* This is a general integral Y(i,j)=∫(ξ^i)(ϕ^j)\[DifferentialD]ξ *)

Y[n_,i_,j_]:=NIntegrate[x^i Evaluate[ϕSol[3][x]]^j,{x,0,ξ1[n]}]

(* For example, Y[2,3] is obtained when calculating M=4π∫(r^2)ρ\[DifferentialD]r *)

Y[3,2,3]

(* Interesting Indices *)

ξ1[1.5]

indices={0,1,1.5,2,3,4};

(* The usual polytropic table *)

MatrixForm[Prepend[Table[{indices[[k]],Evaluate[ξ1[indices[[k]]]],Y[indices[[k]],2,3]},{k,1,Length[indices]}] ,{"n","(Subscript[ξ, 1])","Y[2,3]"}]
]

(* This is a table for n=3, with different i and j *)

ivalues={0,1,2,3,4,5,6};jvalues={0,1,2,3,4,5,6,7,8,9,10,15,20,25};

restemp=Table[Y[ivalues[[i]],jvalues[[j]]],{i,1,Length[ivalues]},{j,1,Length[jvalues]}];

MatrixForm[restemp]

res=Prepend[Transpose[Prepend[restemp,Table[Style["j="<>ToString[jvalues[[j]]],Underlined],{j,1,Length[jvalues]}]]],Table[If[i==0,"",Style["i="<>ToString[ivalues[[i]]],Underlined]],{i,0,Length[ivalues]}]];

MatrixForm[res]

(* Y(i,j)=∫(ξ^i)(ϕ^j)\[DifferentialD]ξ *)

(* Radius of Gyration squared is *)

αG2[n_]:=((2 Y[n,4,n])/(3 ξ1[n]^2 Y[n,2,n]))

αG2[3]