Strongly related MMa.SE thread: "
How to integrate functions of linearly interpolated data?"
In short, you need an
InterpolatingFunction computed with at least 50 digits of
Precision and use the
"InterpolationPointsSubdivision" preprocessor of
NIntegrate with appropriate
"MaxSubregions" suboption:
f1 = NDSolve[{y''[x] + Sin[y[x]] == 0, y[0] == 0, y[10] == 0}, y, x,
Method -> {"Shooting",
"StartingInitialConditions" -> {y[0] == 0, y'[0] == 2}},
WorkingPrecision -> 50][[1, 1, 2]];
NIntegrate[f1'[t]*f1[t], {t, 0, 5}, WorkingPrecision -> 50,
Method -> {"InterpolationPointsSubdivision", "MaxSubregions" -> 1 + Length[First@f1["Coordinates"]]}]
4.7667705747559250585699457377132228611164733033539
