What you are asking for is a common request, and it can be implemented, a little awkwardly. However, what is being asked for in this situation is seldom what people turn out to want.
Every attempted ode45() step involves execution of ode function 6 times with a variety of times and boundary conditions, using 5 of the values to make a prediction and then trying at the predicted location to see if the prediction was good. If the cross-checking of tolerances succeeds then ode45 "accepts" the step and saves the prediction, and then moves on to the next step starting from the prediction. Every once in a while at unspecified times, ode45 uses the recorded values to synthesize a result without having executed at that exact location, and adds the synthesized values to the output.
If the prediction cross-checking fails then ode45 rejects the proposed step, reduces the step size, and tries again with the more conservative step. Several rejections in a row can happen.
Notice that the output values returned to the user might be all synthetic, never having tested that exact location at all. So there might not be Rx, Ry, J corresponding to the locations that were output.
And notice that even in the best case, there are 5 or more "off track" executions whose results would not be recorded. "I know that slightly up the hill is not going to be the actual path, but testing there gives me useful information about how the track is curving"
There are thus a bunch of results to potentially record corresponding to locations used only for analysis. And as described above, there are output locations for which Rx Ry J were never calculated.
You can record all of the values at all of the locations actually executed for, but unless you are trying to do something like calculate confidence intervals, would that be of any use to you?
I suggest that you create a helper function that accepts a column vector of times and an array of boundary conditions, and which uses vectorized calculations to compute the results, and return everything of interest in separate variables.
Now take your dynamics function and have it call the helper, passing it the current time and boundaries — so all the real work is done by the helper.
Now run the ode45 normally. When you finish, take the two output arrays and pass them to the helper function, getting out the Rx, Ry, J for each of the locations that was actually returned. No global variables needed, no loops, just matrix calculations. No values returned for locations used for information, or rejected locations.
