This work presents a numerical method for evaluating affine relaxations of the solutions of parametric ordinary differential equations. This method is derived from a general theory for the construction of a polyhedral outer approximation of the reachable set (“polyhedral bounds”) of a constrained dynamic system subject to uncertain time-varying inputs and initial conditions. This theory is an extension of differential inequality-based comparison theorems. The new affine relaxation method is capable of incorporating information from simultaneously constructed interval bounds as well as other constraints on the states; not only does this improve the quality of the relaxations but it also yields numerical advantages that speed up the computation of the relaxations. Examples demonstrate that tight affine relaxations can be computed efficiently with this method.