  Extreme mass ratio inspirals (EMRIs), where a compact object orbits a massive
black hole, are a key source of gravitational waves for the future Laser
Interferometer Space Antenna (LISA). Due to their small mass ratio, ($\epsilon
\sim 10^{-4}$--$10^{-7}$), the binary evolves slowly and EMRI signals will be
in-band for years. Additionally, astrophysical EMRIs are expected to have
complex dynamics featuring both spin-precession and eccentricity. A standard
approach to modelling these inspirals is via the method of osculating geodesics
(OG) which we employ along with a toy model for the gravitational self-force.
Using this method requires resolving tens of thousands radial and polar orbital
librations over the long duration of the signal which makes the inspiral
trajectory expensive to compute. In this work we accelerate these calculations
by employing Near-Identity (averaging) Transformations. However, this averaging
technique breaks down at orbital resonances where the radial and polar
frequencies are an integer ratio of each other. Thus, we switch to a partial
averaging transformation in the vicinity of the resonance where the dynamics
are characterised by the slow evolution of the so-called "resonant phase".
Additionally, we develop an optimal switching criterion to minimise the
computation time while maximising accuracy. We find the error in the waveform
phase is improved from $\mathcal{O}(\epsilon^{-1/2})$ in the fully averaged
scheme to $\mathcal{O}(\epsilon^{4/7})$ in the switching scheme. At the same
time, this scheme improves the scaling of the computation time from being
inversely proportional to $\epsilon$ using OG, to a very weak scaling with
$\epsilon$. This results in a speed-up of at least two orders of magnitude for
LISA EMRIs with room for further optimisation.
