Discrete gradient methods are a class of numerical integrators producing solutions with exact preservation of first integrals of ordinary differential equations. In this paper, we apply order theory combined with the symmetrized Itoh--Abe discrete gradient and finite differences to construct an integral-preserving fourth-order method that is derivative-free. The numerical scheme is implicit and a convergence result for Newton's iterations is provided, taking into account how the error due to the finite difference approximations affects the convergence rate. Numerical experiments verify the order and show that the derivative-free method is significantly faster than obtaining derivatives by automatic differentiation. Finally, an experiment using topographic data as the potential function of a Hamiltonian oscillator demonstrates how this method allows the simulation of discrete-time dynamics from a Hamiltonian that is a combination of data and analytical expressions.

18 pages, 7 figures