The paper demonstrates symmetric integral operator and interpolation based numerical approximations for linear and nonlinear ordinary differential equations (ODEs) with oscillatory factor x′=ψ(x)+χω(t), where the function χω(t) is an oscillatory forcing term. These equations appear in a variety of computational problems occurring in Fourier analysis, computational harmonic analysis, fluid dynamics, electromagnetics, and quantum mechanics. Classical methods such as Runge–Kutta methods etc. fail to approximate the oscillatory ODEs due the existence of oscillatory term χω(t). Two types of methods are presented to approximate highly oscillatory ODEs. The first method uses radial basis function interpolation, and then quadrature rules are used to evaluate the integral part of the solution equation. The second approach is more generic and can approximate highly oscillatory and nonoscillatory initial value problems. Accordingly, the first-order initial value problem with oscillatory forcing term is transformed into highly oscillatory integral equation. The transformed symmetric oscillatory integral equation is then evaluated numerically by the Levin collocation method. Finally, the nonlinear form of the initial value problems with an oscillatory forcing term is converted into a linear form using Bernoulli’s transformation. The resulting linear oscillatory problem is then computed by the Levin method. Results of the proposed methods are more reliable and accurate than some state-of-the-art methods such as asymptotic method, etc. The improved results are shown in the numerical section.