We present a numerical scheme for an efficient discretization of nonlinear systems of differential equations subjected to highly oscillatory perturbations. This method is superior to standard ODE numerical solvers in the presence of high frequency forcing terms,and is based on asymptotic expansions of the solution in inverse powers of the oscillatory parameter w, featuring modulated Fourier series in the expansion coefficients. Analysis of numerical stability and numerical examples are included.