ClearAll["Global`*"]
Off[NDSolve::mxst];
Off[NDSolve::ndsz];
G = 6.673*10^-8;
\[Rho]c = 1;
X = {0.7, 0.3};
\[Mu]x = {2, 4/3};
\[Mu] = 1/Total[X/\[Mu]x];
mp = 1.67262158*10^-24;
kb = 1.3806503*10^-16;
Pc = (\[Mu] \[Rho]c)/mp kb Tc;
lef[n_?NumericQ] := \[Theta][\[Xi]] /.
   Module[{\[Epsilon] = $MachineEpsilon},
     NDSolve[{1/\[Xi]^2 D[\[Xi]^2 D[\[Theta][\[Xi]], \[Xi]], \[Xi]] + \
\[Theta][\[Xi]]^n == 0, \[Theta][\[Epsilon]] ==
        1, \[Theta]'[\[Epsilon]] ==
        0}, \[Theta][\[Xi]], {\[Xi], \[Epsilon], 30},
      MaxSteps -> 10000]][[1, 1]];
leterm[n_?NumericQ] := \[Xi] /.
   Module[{func = Re[lef[n]], \[Epsilon] = $MachineEpsilon},
    FindRoot[func == 0, {\[Xi], 1, 2}]];
ilef[x_?NumericQ, n_?NumericQ] :=
  Module[{func = lef[n], \[Epsilon] = $MachineEpsilon},
   NIntegrate[func, {\[Xi], \[Epsilon], x}]];
elef[x_?NumericQ, n_?NumericQ] :=
  Module[{func = lef[n], \[Epsilon] = $MachineEpsilon,
    xmax = leterm[n]}, NIntegrate[func, {\[Xi], x, xmax}]];
dr = Sqrt[((n + 1) Pc)/(4 \[Pi] G \[Rho]c^2)];
r = \[Xi] dr;
volilef[x_?NumericQ, nn_?NumericQ, TTc_?NumericQ] :=
  4 \[Pi] Module[{func = lef[nn], \[Epsilon] = $MachineEpsilon},
    NIntegrate[
     r^2 dr func /. {n -> nn, Tc -> TTc}, {\[Xi], \[Epsilon], x}]];
le[\[Xi]\[Xi]_?NumericQ, n_?NumericQ] :=
  lef[n] /. \[Xi] -> \[Xi]\[Xi];