在尝试求解把vk/vkdig这样的变上限积分携带自变量代入到微分方程中求解时，发现它的计算速度非常的慢。请问该怎么优化呢？

完整代码如下：
Mpl = 1;(*Planck单位制*)
g = 106.75;(*SM自由度*)
n = 2;
V0 = 10^-10;
\[Phi]0 = 0.1;
c = -1.5;
p = 0;
CY = 3*10^-9;
\[Rho]r0 =
6.18*10^-14;(*0.1/6.18e-14 0.15/1.17e-13 0.2/1.83e-13*)
endtime = 7*10^6;
V = V0/2 (\[Phi][\[Eta]]/Mpl)^2;
\[ScriptCapitalH] = Sqrt[
a[\[Eta]]^2/(3 Mpl^2) (\[Rho]\[Phi] + \[Rho]r[\[Eta]])];
Y = CY*T^c \[Phi][\[Eta]]^p Mpl^(1 - p - c);
\[Rho]\[Phi] = \[Phi]'[\[Eta]]^2/(a[\[Eta]]^2 2) + V;
(*\[Rho]r[t]=(\[Pi]^2*g*T[t]^4)/30;*)
T = ((30 \[Rho]r[\[Eta]])/(\[Pi]^2*g))^(1/4);
Q = (a[\[Eta]] Y)/(3 \[ScriptCapitalH]);
eom\[Phi] = \[Phi]''[\[Eta]] + (2 \[ScriptCapitalH] +
a[\[Eta]] Y) \[Phi]'[\[Eta]] +
a[\[Eta]]^2 D[V, {\[Phi][\[Eta]], 1}];(*\[Phi]的运动方程*)
eomr = D[\[Rho]r[\[Eta]], \[Eta]] +
4 \[ScriptCapitalH]*\[Rho]r[\[Eta]] - (Y*\[Phi]'[\[Eta]]^2)/
a[\[Eta]];(*Subscript[\[Rho], r]的运动方程*)
eoma = a'[\[Eta]]/a[\[Eta]]^2 - Sqrt[
1/(3 Mpl^2) (\[Rho]\[Phi] + \[Rho]r[\[Eta]])];(*a的运动方程*)
{re\[Phi], \[Rho]reff, aeff} =
NDSolveValue[{eom\[Phi] == 0, eomr == 0, eoma == 0,
a[0] == 1, \[Phi][0] == \[Phi]0, \[Phi]'[0] ==
0, \[Rho]r[0] == \[Rho]r0}, {\[Phi], \[Rho]r, a}, {\[Eta], 0,
endtime}];
function[k_] := Module[{
mx = 10^-7},
mxeff[\[Eta]_] := Sqrt[aeff[\[Eta]]^2 mx^2];(*代入a[\[Eta]]的结果*)
\[Omega]k[\[Eta]_] := Sqrt[k^2 + mxeff[\[Eta]]^2];
vk[\[Eta]_?NumericQ] :=
1/Sqrt[2 \[Omega]k[\[Eta]]]
Exp[-I*NIntegrate[\[Omega]k[t], {t, 0, \[Eta]}]];
vkdig[\[Eta]_?NumericQ] :=
1/Sqrt[2 \[Omega]k[\[Eta]]]
Exp[I*NIntegrate[\[Omega]k[t], {t, 0, \[Eta]}]];
{D[\[Beta]k[\[Eta]], \[Eta]] -
D[\[Omega]k[\[Eta]], \[Eta]]/(2 \[Omega]k[\[Eta]]) vk[\[Eta]]/
vkdig[\[Eta]]*\[Alpha]k[\[Eta]] == 0,
D[\[Alpha]k[\[Eta]], \[Eta]] -
D[\[Omega]k[\[Eta]], \[Eta]]/(2 \[Omega]k[\[Eta]]) vkdig[\[Eta]]/
vk[\[Eta]]*\[Beta]k[\[Eta]] == 0, \[Alpha]k[0] ==
1, \[Beta]k[0] == 0}
]
\[Beta]keff[k_] :=
NDSolveValue[
function[k], {\[Alpha]k, \[Beta]k}, {\[Eta], 0, endtime},
MaxSteps -> Infinity][[2]]
fx[k_, \[Eta]_] := Abs[\[Beta]keff[k][\[Eta]]]^2
s[\[Eta]_] = ParallelTable[fx[10^i, \[Eta]], {i, -5.6, -5.6, 0.1}];

你这aeff的形式还是挺简单的，在NIntegrate里塞Method -> {Automatic, "SymbolicProcessing" -> 0}应该是个比较快当的改法。更快但更高级的方法应该是把InterpolatingFunction里的函数值列表搞出来拿去参与运算最后再重新构造一个函数。参《What's inside InterpolatingFunction[{{1., 4.}}, <>]》（编号28337）。
另外请注意Integrate的不定积分语法是可以直接处理InterpolatingFunction的，这个性质也可以用来简化编程。

不是，我的意思是Integrate的不定积分语法可以直接处理插值函数（没想到吧，加减乘除不行，积分反而行）：
f = Interpolation[Sin@Range[0, 2 Pi, 0.1]];
expr = Integrate[f[t], t]
intf = expr // Head
ListPlot@intf
注意最后一行用了个目前文档里没写的语法：ListPlot可以直接画InterpolatingFunction。