  We present a novel numerical method for solving ODEs while preserving
polynomial first integrals. The method is based on introducing multiple
quadratic auxiliary variables to reformulate the ODE as an equivalent but
higher-dimensional ODE with only quadratic integrals to which the midpoint rule
is applied. The quadratic auxiliary variables can subsequently be eliminated
yielding a midpoint-like method on the original phase space. The resulting
method is shown to be a novel discrete gradient method. Furthermore, the
averaged vector field method can be obtained as a special case of the proposed
method. The method can be extended to higher-order through composition and is
illustrated through a number of numerical examples.
