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, ε ~ 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 O(ε-1/2) in the fully averaged scheme to O(ε4/7) in the switching scheme. At the same time, this scheme improves the scaling of the computation time from being inversely proportional to ε using OG, to a very weak scaling with ε. This results in a speed-up of at least two orders of magnitude for LISA EMRIs with room for further optimisation.