NDSolve will compute a numeric antiderivative.
Example:
ClearAll[f];
f[t_] := Re@Zeta[1/2 + Sqrt[t] I];
Plot[f[t], {t, 1, 1000}]
The antiderivative:
ListLinePlot@NDSolveValue[{y'[t] == f[t], y[1] == 0}, y, {t, 1, 1000}]
Here's the definite integral int[x, y]:
ClearAll[int];
With[{F = NDSolveValue[{y'[t] == f[t], y[1] == 0}, y, {t, 1, 1000}]},
int[x_, y_] := Subtract @@ F[{y, x}]
];
Plot3D[int[x, y], {x, 1, 1000}, {y, 1, 1000}]
For greater accuracy, use
NDSolveValue[{y'[t] == f[t], y[1] == 0}, y, {t, 1, 1000}, InterpolatingOrder -> All]