We present two numerical schemes of high accuracy for one-dimensional dynamical systems. They are modifications of the discrete gradient method and keep its advantages, including stability and conservation of the energy integral. However, their accuracy is higher by several orders of magnitude.