The concern of this paper is in expanding and computing initial-value problems of the form
y
′=
f
(
y
)+
h
ω
(
t
), where the function
h
ω
oscillates rapidly for
ω
≫1. Asymptotic expansions for such equations are well understood in the case of modulated Fourier oscillators
and they can be used as an organizing principle for very accurate and affordable numerical solvers. However, there is no similar theory for more general oscillators, and there are sound reasons to believe that approximations of this kind are unsuitable in that setting. We follow in this paper an alternative route, demonstrating that, for a much more general family of oscillators, e.g. linear combinations of functions of the form e
i
ωg
k
(
t
)
, it is possible to expand
y
(
t
) in a different manner. Each
r
th term in the expansion is
for some
ς
>0 and it can be represented as an
r
-dimensional highly oscillatory integral. Because computation of multivariate highly oscillatory integrals is fairly well understood, this provides a powerful method for an effective discretization of a numerical solution for a large family of highly oscillatory ordinary differential equations.