Evaluate function at "good" points used by Plot in numerical problems

In principle you should be able to do this with FunctionInterpolation, but it's not super well documented.

As an example, the function Sin[x^2] becomes progressively more curved and therefore needs progressively denser sampling:

int = FunctionInterpolation[Sin[x^2], {x, 0, 10}, MaxRecursion -> 20, InterpolationOrder -> 2]

You can inspect the internals of the interpolation function by evaluating:

List @@ int

You'll notice that the x-values are stored in the 3rd argument and the y-values in the 4th. You can plot them like this:

ListPlot[Transpose @ {int[[3, 1]], int[[4, All, 1]]}] 

The sampling density increases with x:

Histogram[int[[3, 1]], 50]

FunctionInterpolation has a number of options that control how it samples the function. You may need to tinker with these to get a good result (especially the MaxRecursion option often needs increasing). You can find them by evaluating:

Options[FunctionInterpolation]

{AccuracyGoal -> Automatic, InterpolationOrder -> 3, InterpolationPoints -> 11, InterpolationPrecision -> Automatic, MaxRecursion -> 6, PrecisionGoal -> Automatic}

This answer provides more detail.

You can integrate a DAE for your function with NDSolve. I used a low-order integration rule to get dense sampling when the second derivative is large in magnitude. I used a low PrecisionGoal so that the number of points would be low, which helps the visualization below show where the sampling is denser. You can change the precision as desired.

approx = NDSolveValue[{
    x'[t] == 1, x[0] == 0,       (* dummy DE *)
    y[t] == Sin[3 t] - Sin[t]},  (* function to integrate *)
    y, {t, 0, 2 Pi}, 
   Method -> {"IDA", "MaxDifferenceOrder" -> 1},
   PrecisionGoal -> 2, AccuracyGoal -> 3];

Plot[Sin[3 t] - Sin[t], {t, 0, 2 Pi},
 Mesh -> {Flatten@approx@"Grid"}, MeshStyle -> Red]

You can get the function values with:

approx@"ValuesOnGrid"

(*  {0.`, 0.0002, ..., -0.0338323, -4.92661*10^-16}  *)