Highly oscillatory ordinary differential equations (ODEs) pose a computational challenge for standard solvers available in scientific computing libraries. These conventional methods are typically based on a polynomial approximation, resulting in there being several timesteps per oscillation period, which leads to runtimes scaling as πͺ(π), with π being the oscillation frequency. This can become prohibitively slow.