A reverse engineering approach is taken to develop efficient numerical techniques for solving nonlinear ODEs through the numerical solution of a particular problem in engineering. The problem is the dynamic analysis of a mechanical system under a severely irregular seismic signal or earthquake load. It is a challenging problem which requires strong formulation to be tackled. Its solution is a rich source of information that could be exploited through artificial intelligence (AI) and statistical computation. This work devises a new strategy to conduct such a task and construct a series of 3∼5 points high-precision time integrators. The integrators belong to single-step two-derivative class of methods. The interconnectivity relationships between the weights of the integrators are formulated by linear regression and heuristics. These relationships make it possible to generate a wide spectrum of ODE-solving techniques in a continuous weighting space. In an innovative attempt, strong Hermite interpolators are constructed to couple with the integrators and improve the accuracy of the presented algorithms. They precisely evaluate the instep points of the integrator. All these issues are appropriately collected in the presented methodology. Four classes of weighting rules are presented for the integrators: 1) Symmetry weighting rule (SWR) 2) Consistency weighting rule (CWR), 3) Fundamental weighting rule (FWR), and 4) Auxiliary weighting rule (AWR). The first two rules are adopted from outward form of the well-configured quadratures available in literature. The third and the most important one, the so-called FWR, is disclosed here for the given integrator. It formulates the correlation between the weights of an optimal integrator. Finding the min-error solution of the vibration equation guides us to the FWR. Evolutionary grey wolf optimizer (GWO), accompanied with a statistical multilinear regression, extracts the FWR. This research pioneers the application of swarm-intelligence in formulation of the FWRs for the given two-derivative integrator. Finally, all the weighting rules and Hermite interpolators are combined into some unified algorithms so-called γβII−q+rP algorithms to efficiently cope with initial value ODEs. Great performance of the presented methods is demonstrated through the numerical examples.