After integrating your ODE-system you will have
for some
, and you have the function for your ODE-equations, f(t,X), one way to calculate the accelerations at times
will be to simply call: dxd2xdyd2ydt(i,:) = f(t_i,X(i,:));
Typically I use a regularly space array for the time in the call to ode45 and use the gradient-function, but the above should work - that is after all what we expect odeNN: to give us a solution X where that is true...
HTH