The most accurate numerical method is to let the integrator decide for the best step sizes. Then you get the optimal solution between the accuracy through small step sizes and reducing the accumulated rounding errors by a small number of steps.
If you need the output of the trajectory to certain time points an interpolation after the integration is a good strategy. I did not take a look into the code of ODE45 for a long time now, but many other professional integrators do not modify the step size of the integration to meet requested time points for the output, but they re-use the calculated trajectory and the intermediate results from the Runge-Kutta steps for an interpolation.
This means, that your point 3. ist the best choice. To increase the accuracy, reduce the tolerances.
Another approach would be define a tspan to obatin the trajectories to defined time points, and to remove the event times afterwards:
% !!!UNTESTED CODE!!!
y0 = 0;
tspan0 = 1:1:10;
t0 = 1;
ready = false;
Opt = odeset('event', ???);
tAll = [];
yAll = [];
while ~ready
tspan = [t0, tspan0(tspan0 > t0)]; % [TYPO fixed]
[t, y] = ode45(@fcn, tspan, y0, Opt);
tAll = cat(1, tAll, t);
yAll = cat(1, yAll, y);
ready = (t(end) == tspan0(end));
if ~ready
y0 = y(end, :);
t0 = t(end);
end
end
Jitter = ~ismembertol(tAll, tspan0);
tAll(Jitter) = [];
yAll(Jitter, :) = [];
So crop tspan in each iteration and remove the irregular time points introduced by the event functions afterwards.
