We investigate positively invariant sets for an ordinary differential equation (ODE), that are also positively invariant for numerical methods to compute its solution. In particular, we show that for an ODE with an exponentially stable equilibrium and an arbitrary compact subset of its basin of attraction, we can establish the existence of a larger compact set that is positively invariant for both the ODE, one-step explicit- and multi-step numerical methods, and even predictor-corrector multi-step methods. We demonstrate in an example the use of this method when computing a contraction metric for an ODE with an exponentially stable equilibrium.