This is my first post on this site, please forgive me if I screw up to copy my code in a nice way..
I am trying to solve a set of 3 coupled differential equations with NDSolve, where one of the terms in the equations is obtained by NIntegrate. The set of equations looks as follows
NDSolve[{D[X[A], A] == (A^(7/2)/TRHb^5)*((Sqrt[3]*Mpl*RF[A,10^8,10^12])/Sqrt[Φ[A] + R[A]/A + (X[A]*Sqrt[mDM^2 + 9*κb^2*(A^(-3/2))^(1/2)])/TRHb]),
D[Φ[A], A] == (-(A^(1/2)/TRHb^2))*((Sqrt[3]*Mpl*ΓbSM)/Sqrt[Φ[A] + R[A]/A + (X[A]*Sqrt[mDM^2 + 9*κb^2*(A^(-3/2))^(1/2)])/TRHb]),
D[R[A], A] == (ΓbSM*Φ[A] - (1/TRHb^4)*(2*Sqrt[mDM^2 + 9*κb^2*(A^(-3/2))^(1/2)]*RF[A,10^8,10^12]*A^3))*((Sqrt[3]*Mpl*A^(3/2))/(TRHb^2*Sqrt[Φ[A] + R[A]/A + (X[A]*Sqrt[mDM^2 + 9*κb^2*(A^(-3/2))^(1/2)])/TRHb])),
X[10^(-10)] == 0,
Φ[10^(-10)] == ΦI,
R[10^(-10)] == 0},
{X, Φ, R}, {A, 10^(-10), 10^4}]}
The function RF is the one which I need to solve with NIntegrate. The final goal is to solve this set of equations for X[A,mDM,λ]. In order to do so, I am trying to solve these equations on a "grid" of different mDM and λ values. Afterwards I fix A and use Interpolation in order to get a nice smooth solution in mDM and λ. But I already struggle to get a solution just for a single combination of mDM and λ (mDM = 10^8, λ = 10^12). I guess the main problem is in the function RF. It is defined as follows
RF[(A_)?NumericQ, (mDM_)?NumericQ, (mt_)?NumericQ, (αG_)?NumericQ, (λ_)?NumericQ]:=
(4*Pi/(2048*Pi^6))*NIntegrate[
Norm[(Sqrt[1-(4*mDM^2)/s]*((mDM^2*s^3*(1-(4*mDM^2)/s))*(4 + 12*αG^2)^2)*κb*(A^(-3/2))^(1/4)*Sqrt[s]*BesselK[1, Sqrt[s]/(κb*(A^(-3/2))^(1/4))])/(λ^4*((s - mt^2)^2 + mt^2*((mt^3*((1/8)*Sqrt[1 - (4*89^2)/mt^2] + (3*αG^2)/16))/(Pi*λ^2) + (mt^3*Sqrt[1 - (4*mDM^2)/mt^2]*mDM^2*(1 - (4*mDM^2)/mt^2))/((Pi*λ^2)*8*mt^2))^2))],
{s, 4*mDM^2, Infinity}, Exclusions -> {mt^2, ma^2}]
As a function of A, RF[A] should look like FIG1. in [ https://arxiv.org/abs/1811.01947 ]. But I fail to reproduce this form in the rage mDM = [10^8 - 10^14]. As soon as NIntegrate hits the point s ~ mt^2, I get problems. I was trying to play around with the NIntegrate settings, but nothing really helped. I obtained the best solutions with "DoubleExponential" and "Trapezoidal", but in both cases I get Singaluarity errors. I guess I need to get first an error free solution in NIntegrate before I can use RF in NDSolve. How can I get a nice integration for the entire range of mDM?
The parameters I used were
ma = 0;αS = 1; αV = 1/Log[10^5]; αH = 1; αG = 1/Log[10^5]; βG = 0; βV = 0; αχ = 1; αχV = 1; βχ = 0; βχ5 = 0; g = 100; V = 10^5; Mpl = 2*10^18; κinf = 10^17; κb = 9*10^8; ΦI = 4*10^43; ΓinfSM = 4*10^8; ΓbSM = 7/10^11; minf = 6*10^14; mt = 2*10^10; TRHinf = 2*10^13; TRHb = 7*10^4;
Update
Here are some plots of the function RF[A] for the parameters mDM = 10^8, 10^9, 10^10, 10^11, 10^12, 10^13. Inside NIntegrate I used
{s, 4 mDM^2, mt^2, Infinity}, Method -> "DoubleExponential", WorkingPrecision -> 10, AccuracyGoal -> 10.
Timing[Plot[RF[A, mDM, mt, αG, λ], {A, 10^(-10), 10^10}, ScalingFunctions -> {"Log", "Log"}, PlotRange -> Automatic, PerformanceGoal -> "Quality"]]
And I get the following error for the case of mDM=10^12
NIntegrate::ncvi: NIntegrate failed to converge to prescribed accuracy after 9 iterated refinements in s in the region {{4.00000000000000000000000000000000000000000000000000000000000*10^20,\[Infinity]}}. NIntegrate obtained 9.2531746437536964660643301435717812896328729569142552371394`60.*^54 and 2.60922345176519777170004482614311324727140650854448718254727`60.*^53 for the integral and error estimates.
I guess before I use RF[A] as an input for NDSolve I need to get it as a smooth function for the entire parameter space.