Perturbed harmonic oscillators describe many models in physics and engineering. A method for solving this
type of problem is based on the utilization of Scheifele’s functions, consisting of a refinement of the Taylor series method.
One disadvantage of the method is that it is difficult to determine, in each case, the recurrence relations necessary for arriving
to the solutions. In this paper we construct the numerical method especially adapted to the integration of oscillators,
using developments in the form of G-functions. In addition, two computer applications related to highly oscillatory problems
are implemented, and recurrence relations are determined in each one. The results show better precision in the application
of G-functions, compared to other known methods implemented in MAPLE V.