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

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} ]


• FYI: your Friedman equation is (most likely) wrong. To quickly see this, evaluate H[0] which doesn't give H0 as expected (ie the present day value of the Hubble parameter). Commented Dec 22, 2023 at 18:18
• Actually I have not derived that .. It is from Baumann's book (eqn-3.169) Commented Dec 22, 2023 at 18:35
• Ok, it's probably an expansion around the matter-radiation equality. Still it might be prudent to use a correct expression and check how the results changes Commented Dec 22, 2023 at 18:54
• Sure thanks.. can you provide some source for the expression?.. Commented Dec 23, 2023 at 5:36
• See the last equation in the section "Density Parameter" here: en.wikipedia.org/wiki/Friedmann_equations You can ignore curvature and use the equality redshift to relate the radiation to matter. I hope this helps. Commented Dec 23, 2023 at 10:55

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]


• Thanks, Can you tell me one more thing? If I want to find the value of z where Xe is 0.5.. How to do that? Commented Dec 22, 2023 at 18:30
• To find $z$ that solves $X_e(z) =.5$, evaluate FindInstance[{s[z] == .5, 10 < z < 1800}, z]. to get (* {{z -> 1376.57}} *)`. Commented Dec 22, 2023 at 19:18
• Thank you once again for your help. Commented Dec 29, 2023 at 1:00