In Mathematica, using WL, your code
s = NDSolve[{y'[x] == y[x] Cos[x + y[x]], y[0] == 1}, y, {x, 0, 30}]
returns a set of solutions; in this case, there is only one solution. First[s] yields that solution. One can plot the solution with
Plot[y[x] /. First[s], {x, 0, 30}]
One can get a single value of y at a specified x, say x=4.5, thus:
y[4.5] /. First[s]
One can get a table of values like this:
Table[y[x] /. First[s], {x, 0, 30, 0.5}]
Note that for values of x between the steps, an interpolation error is added to the integration error, and it usually has a lower order of precision than the integration method. For plotting and other applications, it is usually small enough not to matter. But if you need greater accuracy, add the option InterpolationOrder -> All to the NDSolve[] call.
If you want the actual steps taken by NDSolve, the following extracts the x and y values (separately, which is how they're stored):
xgrid = y["Grid"] /. First[s] (* xgrid = {{x0}, {x1}, {x2}, ...} *)
yvals = y["ValuesOnGrid"] /. First[s] (* yvals = {y0, y1, y2, ...} *)
Use Flatten[xgrid] if you want a flat list {x0, x1, x2,...} of x coordinates. One could also use xgrid = First[y["Coordinates"] /. First[s]] instead of "Grid" and Flatten. The following is are quick ways to plot the steps:
ListPlot[y /. First[s]] (* Not y[x]! *)
ListLinePlot[y /. First[s]] (* Connects the dots between the steps *)
I hope that helps. I don't use the framework shown your question. I don't know how graphics work in it, for instance. I'm just hoping that all you need are the WL codes for evaluating and plotting the solution.