Abstract It is well known that the application of one-step or linear multistep methods to an ordinary differential equation with first integrals will destroy the conserving quantities. With the use of stabilization techniques similar to Ascher, Chin, Reich (Numer. Math. 67 (1997) 131–149) we derive stabilized variants of our original problem. We show that variable step size one-step and linear multistep methods applied to the stabilized equation will reproduce that phase portrait correctly. In particular, this technique will conserve first integrals over an infinite time interval within the local error of the used method.