In this paper, we develop a class of robust numerical methods for solving dynamical systems with multiple time scales. We first represent the solution of a multiscale dynamical system as a transformation of a slowly varying solution. Then, under the scale separation assumption, we provide a systematic way to construct the transformation map and derive the dynamic equation for the slowly varying solution. We also provide the convergence analysis of the proposed method. Finally, we present several numerical examples, including ODE system with three and four separated time scales to demonstrate the accuracy and efficiency of the proposed method. Numerical results verify that our method is robust in solving ODE systems with multiple time scale, where the time step does not depend on the multiscale parameters
