For the integration of the classical equations of motion in the Born-Oppenheimer approach, each time the energy and gradient of the potential energy surface are needed, a properly converged wave function is calculated. If Hessians (second derivatives) can be calculated, significantly larger steps can be taken in the numerical integration of the equations of motion without loss of accuracy. Even larger steps can be taken with a Hessian-based predictor-corrector algorithm. Since updated Hessians are used successfully in quasi-Newton methods for geometry optimization, it should be possible to improve the performance of trajectory calculations using updated Hessians. The Murtagh-Sargent (MS) update, the Powell-symmetric-Broyden (PSB) update and Bofill's update (a weighted combination of MS and PSB) were tested, and Bofill's update was found to be the best. Slightly smaller step sizes were needed with Hessian updating to maintain good conservation of the energy, but this was more than compensated by the reduction in total computational cost. An overall factor of 3 in speed-up was obtained for trajectories of systems containing 4 to 6 heavy atoms computed at the HF/3-21G level.