  We show how to increase the order of one-dimensional discrete gradient
numerical integrator without losing its advantages, such as exceptional
stability, exact conservation of the energy integral and exact preservation of
the trajectories in the phase space. The accuracy of our integrators is higher
by several orders of magnitude as compared with the standard discrete gradient
scheme (modified midpoint rule) and, what is more, our schemes have very high
accuracy even for large time steps.
