# Performing a FindRoot from Numerical integration

I am trying to evaluate:

FindRoot[Xfnew[γ, w, Hi, mx]*92*10^23*γ^(4/(1 + w)) == 1198/10000, {Hi, hs},
MaxIterations -> 10000]


for some starting value hs, and γ and w are known numbers listed below in the definitions. I have to evaluate this FindRoot in a range of mx, and plot it, something like:

LogLogPlot[Hi /. FindRoot[Xfnew[γ, w, Hi, mx]*9.2*10^24*γ^(4/(1 + w)) ==
0.1198, {Hi, 1.6946*10^-8}], {mx, 10^(-4), 10^(-10)}, PlotStyle -> Green]


I have successfully evaluated this in the range of mx (0.1 to 0.001), and (10^-10 to 10^-19), with starting values 0.0046 and 0.1 respectively.

But in the range 10^(-4) to 10^(-10), I am not able to. I tried using a "for loop" and break this range to form a list {mx,Hi} with small steps of mx(and doing a ListLogLogPlot afterwards), and perturbing the initial starting of FindRoot, like this:

For[i = 1, i < 2000, i++,h = hs; h[i] = Hi /. Check[FindRoot[
Xfnew[γ, w, Hi, mxs - i*step]*92.*10^23*γ^(4/(
1 + w)) == 1198./10000, {Hi, h[i - 1]}, MaxIterations -> 10000],
{Hi -> 25}]; If [h[i] == 25, Break[], AppendTo[hin, {mxs - i*step, h[i]}]]]


where hs and mxs are some starting values of Hi and mx, and I break the loop if FindRoot generates error messages. And "step" is some stepsize I have chosen.

But, currently my step size for this loop has decreased so much, that it would take me ages to cover the entire range.I decrease my stepsize, whenever I get the error:

Machine precision is insufficient to achieve the accuracy

Is there a faster way to cover this range? I have tried increasing Working precision, it still doesnot help.

Function definitions

Xfnew is a result of a complicated list of functions listed below, and the values of gamma and w are listed below:: (Xfnew and the values of gamma and w are listed at the End) (I apologize for the Greek letters)

g = 61;
k2 = 25/100;
k1 = 2/10;
ai = 1;
Trh[γ_?NumericQ, Hi_?NumericQ] := k2*γ*(( Hi)^(1/2))
arh[γ_?NumericQ, w_?NumericQ] := ai*γ^(-4/(3 (1 + w)))
Hrh[a_?NumericQ, γ_?NumericQ, w_?NumericQ, Hi_?NumericQ] :=
Hi (a/ai)^(-3 (1 + w)/2)
H[a_?NumericQ, γ_?NumericQ, w_?NumericQ, Hi_?NumericQ] :=
Hi γ^2 (a/arh[γ, w])^-2
Trh11[a_?NumericQ, γ_?NumericQ, w_?NumericQ,
Hi_?NumericQ] := ((k1 (γ Hi)^(1/2))/((1 + 3/5 w)^((1/
4)))) ((a^(-3 (1 - w)/2) - a^-4)^((1/4)))
T[a_?NumericQ, γ_?NumericQ, w_?NumericQ, Hi_?NumericQ] :=
Trh[γ, Hi] arh[γ, w]/a
nXeq[a_?NumericQ, γ_?NumericQ, w_?NumericQ, Hi_?NumericQ,
mx_?NumericQ, T_Symbol] := g/(2 π^2)
mx^2 T[a, γ, w, Hi] BesselK[2, mx/T[a, γ, w, Hi]]
σv0[a_?NumericQ, γ_?NumericQ, w_?NumericQ, Hi_?NumericQ,
mx_?NumericQ, T_Symbol] := 3 π ((mx^2)/
8) (3 BesselK[1, mx/T[a, γ, w, Hi]]^2/(5 BesselK[2, mx/T[a, γ, w, Hi]]^2) + 2/5 +
4/5 T[a, γ, w, Hi]/mx BesselK[1, mx/T[a, γ, w, Hi]]/
BesselK[2, mx/T[a, γ, w, Hi]] + 8/5 T[a, γ, w, Hi]^2/mx^2)
σv1[a_?NumericQ, γ_?NumericQ, w_?NumericQ, Hi_?NumericQ,
mx_?NumericQ, T_Symbol] := π (
mx^2) (11/20 BesselK[1, mx/T[a, γ, w, Hi]]^2/
BesselK[2, mx/T[a, γ, w, Hi]]^2 + 9/20 + 13/20 T[a, γ, w, Hi]/
mx BesselK[1, mx/T[a, γ, w, Hi]]/BesselK[2, mx/T[a, γ, w, Hi]] +
13/10 T[a, γ, w, Hi]^2/mx^2)
σvtotal[a_?NumericQ, γ_?NumericQ, w_?NumericQ,
Hi_?NumericQ, mx_?NumericQ, T_Symbol] :=
4 σv0[a, γ, w, Hi, mx, T] + (45 + 12) σv1[a, γ, w, Hi, mx, T]
g1[a_?NumericQ, γ_?NumericQ, w_?NumericQ, Hi_?NumericQ,
mx_?NumericQ] := (mx*a^2)/(Hrh[a, γ, w, Hi]*Trh[γ, Hi]^3) σvtotal[
a, γ, w, Hi, mx, Trh11] nXeq[a, γ, w, Hi, mx, Trh11]^2;
g2[a_?NumericQ, γ_?NumericQ, w_?NumericQ, Hi_?NumericQ,
mx_?NumericQ] := (mx*a^2)/(H[a, γ, w, Hi]*Trh[γ, Hi]^3) σvtotal[
a, γ, w, Hi, mx, T] nXeq[a, γ, w, Hi, mx, T]^2;
Xfnew[γ_?NumericQ, w_?NumericQ, Hi_?NumericQ, mx_?NumericQ] :=
NIntegrate[g1[a, γ, w, Hi, mx], {a, 1, arh[γ, w]}] +
NIntegrate[g2[a, γ, w, Hi, mx], {a, arh[γ, w], ∞}]
γ = 1/10;
w = 5/10;


Experimentation indicates that a larger value of WorkingPrecision indeed is necessary. However, too large a value is computationally expensive and also sometimes causes the computation to crash. WorkingPrecision -> 30 seems to be a good compromise. In addition, computations for mx in the range of {10^-6, 10^-4} do not converge well without very good initial guesses. An initial guess as a function of mx fit to values of Hi for mx in the vicinity of 10^-2 and of 10^-7 and again a bit of experimentation suggests 0.2 mx^1.7 for mx > 10^-6. Then,
ParallelTable[{mx = N[10^i], 