# How to speed up the production of this plot?

This is a followup question to this question. From these two answers I can build the function sol[r][[1]][t0] as follows:

ClearAll["Global*"]
Md = 10^(-9);
P = 10;
R = 10^4;
α = 10^(-2);
ϵ = 10^(-4);
γ = 10^(-2);
ke = 0.02*(1 + 0.6625);
k0 = 5*10^20;
σ = 5.67*10^-8;
Rg = 8315;
c = 3*10^8;
G = 6.67/10^11;
M = 2.8*10^30;

Ωk[r_] := Sqrt[(G*M)/r^3];
μ = Md/(3*Pi);
κ = (27*ke)/(2*σ) Rg/μ;
ic = {Co[0] == 1, β[0] ==
0, Σ[
0] == (μ^(3/5)*Ωk[
r]^(2/5)) (κ^(-1/5)*α^(-4/5)*Co[0]^(-1/5)),
h[0] == (κ*α*Σ[0]^2*
Co[0]/Ωk[r]^5)^(1/6),
T[0] == (1/2)*(Ωk[r]*h[0])^2*(μ/
Rg)*(1/(1 + β[0])),
Kkr[0] == (k0*(Σ[0]*h[0]))/T[0]^(7/2)};

eq = {Σ'[
t] == -Σ[
t] + κ^(-1/5) α^(-4/5) μ^(3/
5) Ωk[r]^(2/5)*Co[t]^(-1/5),
h'[t] == -h[
t] + (κ α Σ[t]^2 Ωk[
r]^(-5) Co[t])^(1/6),
T'[t] == -T[t] +
1/2 μ/Rg (Ωk[r]^2 h[t]^2)/(1 + β[t]),
Kkr'[t] == -Kkr[t] + k0 Σ[t]/h[t]*T[t]^(-7/2),
β'[
t] == -β[t] + μ/
Rg (4 σ)/(3 c) T[t]^3/Σ[t] h[t],
Co'[t] == -Co[t] + (1 + β[t])^4*(1 + Kkr[t]/ke)};

t0 = 20; sol =
ParametricNDSolveValue[{eq, ic}, {Σ, h,
T}, {t, 0, t0}, {r}]


Now I need to perform the following integral for a plot:

int[r_?NumericQ, t_?NumericQ] := sol[r][[1]][t0] r^3 Cos[2 CapitalOmega]k[r] t]
F1[t_?NumericQ] := NIntegrate[int[r, t], {r, 10^6, 10^8}]
Plot[F1[t], {t, 0, 3*10^3}]


In a whole night this code has not produced the plot. I tried with memoization for the function int:

 int[r_?NumericQ, t_?NumericQ] := int[r,t]=sol[r][[1]][t0] r^3 Cos[2 CapitalOmega]k[r] t]


but nothing has changed. I tried with with ParallelTable, but Mathematica says it cannot communicate with the kernels (or something like that) and a simple Table[F1[t],{t,0,1000,100}] still has not produced an output.

How can I produce the desired plot?

EDIT Using memoization on sol[r][[1]][t0] (h[r_?NumericQ]:=h[r]=sol[r][[1]][t0]) makes things faster, but it still takes about half an hour on my laptop. Is it possible to do better? Moreover there are some numerical issues in the plot:

EDIT2 I have solved the numerical issues by using MinRecursion->5, but it slows down the computation even more. Any suggestions are welcome.

EDIT3 Would it be possible to integrate using NDSolve like in this question? I tried with:

ParametricNDSolveValue[F[t]==sol[r][[1]][600]r^3Cos[ CapitalOmegak[r] t],{r,10^6,10^8},{t}]


but it does not work.

• You have constants that differ by dozens of orders of magnitude, which will make numerical calculations slow and unreliable. Please, work in natural units, where all constants are of order 1. Thanks. – AccidentalFourierTransform Aug 22 '19 at 12:07
• @AccidentalFourierTransform It's not necessary when options about precision and accuracy are proper enough AFAIK, an example: mathematica.stackexchange.com/q/23696/1871 (Compare Goofy and Michael's answers there. ) – xzczd Aug 22 '19 at 12:23
• @xzczd: That's true to a certain extend, but properly re-scaling the problem will make it possible to stay within the realm of machine precision arithmetic, which is generally faster than arbitrary-precision calculations. It may or may not be the magic bullet for this particular problem, but rescaling is almost always a good idea nonetheless. – Sjoerd Smit Aug 22 '19 at 15:31
• @SjoerdSmit Yeah, WorkingPrecision will usually make the code slower, but Mathematica won't switch to arbitrary-precision calculation if we only adjust PrecisionGoal and AccuracyGoal`. (Just check the timing of Goofy and Michael's solutions. ) – xzczd Aug 22 '19 at 15:53