One idea is to integrate once to get the sample points then compute the remaining integrals as sums:
sample = Transpose@
SortBy[First@Last@Reap[
NIntegrate[x (c = Cos[10 x + Cos[x]]), {x, 0, 5},
EvaluationMonitor :> Sow[{x, c}]]], #[[1]] &];
wt = ((#[[3]] - #[[1]])/2) & /@ Partition[Join[{0}, sample[[1]], {5}], 3, 1];
wt.(sample[[2]] #) & /@ {sample[[1]], sample[[1]]^2, sample[[1]]^3}
{0.0133333, 0.133275, 0.861541}
Note the accuracy is not terribly good, NIntegrate gives:
{0.0125266, 0.131514, 0.855716}
Somethings a bit off in my quick&dirty trapezoid integration but i think this can be made to work.
roll your own
For this example there really is little benefit to NIntegrate's adaptive sampling so we might as well just use a uniform sampling:
np = 651;(*assumed odd for simpsons rule*)
a=5
b=0
wt = ( (a-b)/(np - 1) )/3 Join[{1}, Flatten@ConstantArray[{4, 2}, (np - 1)/2 - 1], {4, 1}] // N;
x = b + (Range[0, np - 1] (a-b)/(np - 1)) // N;
fast = Cos[10 # + Cos[#]] & /@ x // N;
(# fast).wt & /@ {x, x^2, x^3}
{0.0125266, 0.131514, 0.855716}