NIntegrate a vector function without using Indexed and NDSolveValue?

Question

I encoutered an issue about integrand being non-numerical while trying to integrate a vector function using NIntegrate as shown below:

ClearAll[f, y, dummy1, dummy2];
f[s_?NumericQ] := {Sin[s], Cos[4*s]};
NIntegrate[f[s], {s, 0, 3}]

NIntegrate::inum: Integrand f[s] is not numerical at {s} = {0.023872}.

Of course, if I was able to directly plug in the vector itself, everything is good (this is in my scenario impossible since f is actually based on purely numerical results from NIntegrate and NDSolveValue - so I can only use SetDelayed).

I successfully used Indexed as a workaround and also NDSolveValue as mentioned in this answer and both actually solve the problem:

NIntegrate[{Sin[s], Cos[4*s]}, {s, 0, 3}] // AbsoluteTiming
(* {0.00767999, {1.98999, -0.134143}} *)
dummy1[s_?NumericQ] := Indexed[f[s], 1];
dummy2[s_?NumericQ] := Indexed[f[s], 2];
NIntegrate[{dummy1[s], dummy2[s]}, {s, 0, 3}] // AbsoluteTiming
(* {0.00453236, {1.98999, -0.134143}} *)
NDSolveValue[{y'[s] == f[s], y[0] == {0, 0}}, y[3], {s, 0, 3}] // AbsoluteTiming
(* {0.00260897, {1.98999, -0.134143}} *)

Surely, I do have a working solution now. But I wonder if it is somehow still possible to just use NIntegrate with adjusted options in order to avoid Indexed (may become annoying for high dimensions) and NDSolveValue (I noticed for some dummy integrands f[s], there was already a factor of 5 in runtime between the simple NIntegrate and NDSolveValue approach)?

UPDATE: Another reason why I possibly prefer NIntegrate over NDSolveValue is that I may potentially use something like this which might speed things up a lot.