Many scientific and engineering applications require the solution of large systems of initial value problems arising from method of lines discretization of partial differential equations. For systems with widely varying time scales, or with complex physical dynamics, implicit time integration schemes are preferred due to their superior stability properties. These schemes solve at each step linear systems with matrices formed using the Jacobian of the right hand side function. For large applications iterative linear algebra methods, which make use of Jacobian-vector products, are employed. This paper studies the impact that the method of computing Jacobian-vector products has on the overall performance and accuracy of the time integration process. The analysis shows that the most beneficial approach is the direct computation of exact Jacobian-vector products in the context of matrix-free time integrators. This approach does not suffer from approximation errors, reuses the parallelism and data distribution already present in the right-hand side vector computations, and avoids storing or operating on the entire Jacobian matrix.