You could generate a Table of points using the interpolating function;
however, it would be difficult to sample the function at the best
values. If you let Mathematica Plot the function then you can leverage
its sampling algorithm.
Clear[nDSolve]
nDSolve[eqn_, depVar_Symbol, iter_List, opts_: {}] :=
Module[{sol, plt, pts},
sol = depVar /. NDSolve[eqn, depVar, iter,
FilterRules[opts, Options[NDSolve]]][[1]];
plt = Plot[sol[iter[[1]]], iter,
Evaluate[FilterRules[opts, Options[Plot]]]];
pts = Cases[plt, Line[pts_] :> pts, Infinity][[1]];
{sol, plt, pts}]
sol = nDSolve[{y'[x] == y[x] Cos[x + y[x]], y[0] == 1},
y, {x, 0, 30}, PlotRange -> All];
The result includes the interpolating function
Plot[sol[[1]][x], {x, 0, 30}, PlotRange -> All]
the plot used to capture the points
Show[sol[[2]]]
and the list of points
ListLinePlot[sol[[3]], PlotRange -> All]
Bob Hanlon