We introduce a method for eliminating the truncation error produced when numerically integrating an initial value problem using Runge-Kutta-based algorithms. We propose a methodology for constructing an optimal state-space representation that gives zero local numerical truncation error, and in this sense, is the optimal state-space representation for modeling given phase-space dynamics. To that end, we utilize a simple transformation of the state-space equations into their variational form. This process introduces an inherent freedom, similar to the gauge freedom in physics. We then utilize the gauge function to eliminate the numerical truncation error. We show that by choosing an appropriate gauge function the numerical integration error dramatically decreases and one can achieve much better accuracy compared to the standard state variables for a given time-step. Moreover, we derive general expressions yielding the optimal gauge functions given a Newtonian one degree-of-freedom ODE. For the n degrees-of-freedom case we describe a MATLAB® code capable of finding the optimal gauge functions and integrating the given system using the gauge-optimized integration algorithm. In all of our illustrating examples, the gauge-optimized integration outperforms the integration using standard state variables by a few orders of magnitude.