This problem goes back at least to V11. See https://mathematica.stackexchange.com/a/126361 for a couple of workarounds.
You may be thinking that the following implies NIntegrate could handle vector-valued function:
NIntegrate[{Cos[x^3], Sin[x^3]}, {x, 0, 10}]
(* {0.776098, 0.444613} *)
But the following shows that NIntegrate[] is threaded over the list of integrands, and each integrand is integrated separately.
Trace[
NIntegrate[{Cos[x^3], Sin[x^3]}, {x, 0, 10}],
TraceDepth -> 1]
(*
{NIntegrate[{Cos[x^3],Sin[x^3]},{x,0,10}],
{NIntegrate[Cos[x^3],{x,0,10}], NIntegrate[Sin[x^3],{x,0,10}]},
{0.776098,0.444613}}
*)
BTW, @Gianluca's Integrate in effect integrates the interpolating formulas, more or less exactly up to machine precision round off.
