I would program this using an adaptive step size approach where the step size decreases as some power (square root works well, empirically) of the error associated with an assumed period. Specifically:
Start with a point $y_0$ in the dynamics space that is "near" the asymptotic periodic orbit, some time step size $\Delta_0$ which need not be as small as would be required for your eventual desired precisions, and some estimate $\tau_0$ of the period. Divide that period by the step size to find $n_1 = \left\lceil \frac{\tau_0}{\Delta_0} \right\rceil$, and simulate using $n_1$ steps of size $\frac{\tau_0}{n_1}$.
At this point the simulated trajectory will take you to $y_1$ and let the previous state for one less and subsequent state for one more time step be $x_1$ and $z_1$, respectively. Now let $y_1$ be the point of closest approach of the line $x_1z_1$ to the point $x_0$, and let that point be at $\tilde{y}_1 = \xi_1 x_1 + (1-\xi_1 z_1)$. Then your next guess for the period will be
$$t_1 = t_0 + (2\xi_1 - 1) \Delta_0$$
and your degree of regard for the solution is $r_0 = |y_1-y_0|$ for some suitably chosen norm function.
Repeat this step, but every $k$-th iteration (for some suitably chosen $k$, balancing the cost of changing to a different time step size in the dynamic simulation against the desire to adjust the step size aggressively) modify your step size, doing
$$
\Delta_{nk} = \Delta_{(n-1)k} \sqrt{\frac{r_{nk}}{r_{(n-1)k}}}
$$
As I said, the square root might not be the best choice for how aggressively to decrease the step size; for well behaved systems one might do some theoretical work on what the best exponent should be, but $frac12$ seems like a good heuristic choice.
This scheme has the advantage that you can allow the trajectory to start some small distance away from the periodic orbit, regarding the imperfection in early periodicity to be comparable to the errors already inherent in the iterative approximations to a perfect time period. It also, once you have gotten close to the actual orbit and period, converges geometrically, which is much better than just going around for multiple periods.
If the Matlab routine does not let you recover the state for a step before and a step after the specified time interval, you may have to stop one interval early and move forward the last two intervals by explicit linear approximation to the equations.
