Solving a ODE (Peeble's Equation) with initial condition a the maximum input range

Question

I have a differential equation which is, $$ \frac{d X_e}{dz} = \frac{C_r(z)}{H(z)}\times \frac{1}{(1+z)}\Big[\frac{38\eta T^3}{25\pi^2}\alpha (T)X_e^2 - b(T)(1 - X_e) \Big] $$ Where $$ \alpha(T) = \langle \sigma v \rangle \approx 9.8 \Big( \frac{\alpha_s}{m_e} \Big)^2 \Big( \frac{E_I}{T} \Big)^{1/2} \ln\Big( \frac{E_I}{T} \Big)$$

$$ b(T) = \alpha(T) \Big( \frac{m_e T}{2\pi} \Big)^{3/2} e^{-E_I/T}$$ $$C_r = \frac{\Gamma_{2s} + 3\mathcal{P}\Gamma_{2p}}{\Gamma_{2s} +3\mathcal{P}\Gamma_{2p} + \beta }$$ $$\beta = b(T)e^{3E_I/4T}$$

with $\Gamma_{2s} = 8.227 s^{-1}$ and $3\mathcal{P}\Gamma_{2p} = \frac{675}{2432}\Big(\frac{E_I}{T}\Big)^3 \frac{H(T)}{(1 - X_e)\zeta(3)}$. Also, $$ H(z) = \sqrt{\Omega_m} H_0 (1 + z)^{3/2} \Big( 1 + \frac{1 + z}{1 + z_{eq}} \Big) $$ $\Omega_m = 0.309$, $H_0 = 1.5\times 10^{-33}$ and $z_{eq} = 3395$. Finally $T = 0.2329meV\times (1 + z)$

This is peeble's equation in Recombination. I want to solve this in between $z= 10$ to $z= 1800$ with the initial condition $X_e[1800] = 1$. But I am unable to do show. The is an error (see at the end). The code is:

Subscript[E, I] = 13.6;
Subscript[m, e] = 5.109989461*10^5;
\[Eta] = 6*10^-10;
Subscript[z, eq] = 3395;
Subscript[H, 0 ] = 1.5*10^-33;
prefac\[Alpha] = (9.8/(137*Subscript[m, e])^2);
prefac\[Beta] = (Subscript[m, e]/(2*Pi))^(3/2);
Subscript[\[CapitalOmega], m] = 0.3089;

All constants are here.

T[z_] := 23.29*(1 + z)*10^-5;
\[Alpha][z_] := prefac\[Alpha]*Sqrt[Subscript[E, I]/T[z]]*Log[Subscript[E, I]/T[z]];
b[z_] := \[Alpha][z]* prefac\[Beta]*T[z]^(3/2)*E^(-Subscript[E, I]/T[z]);
H[z_] := Subscript[H, 0 ]*Sqrt[Subscript[\[CapitalOmega], m]]*(1 + z)^(3/2)*(1 + ((1 + z)/(1 + Subscript[z, eq])))
Subscript[\[CapitalGamma], s] = 8.227;
Subscript[\[CapitalGamma], p][z_] := (675/2432)*(Subscript[E, I]/T[z])^3*H[z]/((1 -Subscript[X, e][z])*Zeta[3]);
\[Beta][z_] := b[z]*E^(3 Subscript[E, I]/(4 T[z]));
Subscript[C, r][z_] := (Subscript[\[CapitalGamma], s] + Subscript[\[CapitalGamma], p][z])/(Subscript[\[CapitalGamma],s] + Subscript[\[CapitalGamma], p][z] + \[Beta][z]);
s = NDSolve[{Subscript[X, e]'[x] == (1/(H[x]*(1 + x)))*((38*\[Eta]*(T[x])^3*\[Alpha][x]*(Subscript[X, e][x])^2)/(23*Pi^2) - b[x]*(1 - Subscript[X, e][x])), Subscript[X, e][1800] == 1}, Subscript[X, e], {x, 10, 1800} ]

Can someone please help?. If the code is hard to read, here is an image.

Answer 1

The error message means MMA is trying to find the derivative at the initial condition, $x=1800$. Take a look at H[1800] or H[x] /. x->1800: $$\text{6.371926864638378$\grave{ }$*${}^{\wedge}$-29} \left(\frac{1801}{1800_{\text{eq}}+1}+1\right)$$

The problem is that MMA doesn't put in $z_{eq}$. The expert MMA programmers advise us to avoid the use of subscripts. For example see Item 3 of these pitfalls. Note that you have defined the $E_I = 13.6$. That looks like Euler's number with the imaginary unit as a subscript to MMA. That definition won't go away with Clear. See this answer for how to unset subscripts.

Your code runs fine without subscripts:

ClearAll["Global`*"]

EI = 13.6;
me = 5.109989461*10^5;
η = 6*10^-10;
zeq = 3395;
H0 = 1.5*10^-33;
prefacα = (9.8/(137*me)^2);
prefacβ = (me/(2*Pi))^(3/2);
Ωm = 0.3089;

T[z_] := 23.29*(1 + z)*10^-5;
α[z_] := prefacα*Sqrt[EI/T[z]]*Log[EI/T[z]];
b[z_] := α[z]*prefacβ*T[z]^(3/2)*E^(-EI/T[z]);
H[z_] := H0*Sqrt[Ωm]*(1 + z)^(3/2)*(1 + ((1 + z)/(1 + zeq)))
Γs = 8.227;
Γp[z_] := (675/2432)*(EI/T[z])^3*H[z]/((1 - Xe[z])*Zeta[3]);
β[z_] := b[z]*E^(3 EI/(4 T[z]));
Cr[z_] := (Γs + Γp[z])/(Γs + Γp[z] + β[z]);

ode = Xe'[x] == (1/(H[x]*(1 + x)))*
       ((38*η*(T[x])^3*α[x]*(Xe[x])^2)/(23*Pi^2) -
       b[x]*(1 - Xe[x]));

ode /. x -> 1800  /. Xe[1800] -> 1

(* Xe'[1800] == 1.67186 *)

And the solution is

s = NDSolveValue[{ode, Xe[1800] == 1}, Xe, {x, 10, 1800}];
Plot[s[x], {x, 10, 1800}, GridLines -> Automatic]