Solving ODEs numerically while preserving a first integral