Numerical integration of ODEs while preserving all polynomial first integrals

We present a novel method for solving ordinary differential equations (ODEs) while preserving all polynomial first integrals. The method is essentially a symplectic Runge-Kutta method applied to a reformulated version of the ODE under study and is illustrated through a number of examples including Hamiltonian ODEs, a Nambu system and the Toda Lattice. When applied to certain Hamiltonian ODEs, the proposed method yields the averaged vector field method as a special case.