How to speed up the production of this plot?

Question

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.