We present a class of numerical schemes that are derived from a perturbated interpolant that has a varying order polynomial as the assumed solution to the Ordinary Differential Equation (ODE) models on climate change. These numerical integrators are capable of solving problems arising from chemical kinetics, population models, mechanical oscillations, planetary motions, electrical networks, nuclear reactor control, tunnel switching problems, reversible enzyme kinetics. But in this work, we desire to apply them to the numerical solutions of Autonomous Initial Value Problems (IVPs) in ODEs on climate change. At the end we conduct a numerical experiment, the resulting methods, algorithms, and solutions will be predictive tools for the study of climatic change models with applications to big data.