Many important dynamical systems can be modeled one-dimensional non-linear oscillators [1,2]. For information on the properties of such systems, a variety of analytical methods exist which can be used to construct approximations to the periodic solutions [3,4]. However, when very accurate solutions are needed, the usual practical procedure is to numerically integrate the equations of motion (EOM) [5,6]. A major di$culty with numerical techniques is that they can give rise to numerical instabilities; these are solutions of the discrete equations, used for the numerical integration, that do not correspond to any solution of the di!erential equation [5,6]. For conservative oscillators, numerical instabilities arise when the second order "nite-di!erence schemes do not possess a discrete version of the "rst integral that exists for the EOM [7,8]. This "rst integral, for conservative systems, is the energy function [2].
In an earlier paper [8], it was shown how to construct a discrete energy function such that the derived EOM provided an accurate discrete model of the original di!erential equation. However, this paper [8] only derived the EOM for a one-dimensional oscillator for which the potential energy was a quartic function of the dependent variable x, i.e.,