A method of numerical integration is described which was developed for the computation of ballistic missile trajectories. But the technique is generally applicable to the initial-value problem for systems of ordinary differential equations. The associated variational equations are integrated simultaneously to provide the sensitivity coefficients of the solution. This facilitates the inclusion of an extra derivative in the numerical integration. The differences and advantages of the method as compared to standard integration techniques are discussed. A Newton-Raphson iterative solution of a single-step fourth-order Hermite corrector provides unconditional numerical stability. This very favorable property should make the method useful for stiff systems of ordinary differential equations. Polynomials are derived for interpolating the solution and its sensitivity coefficients, and also for extrapolating the solution to provide a predictor for the Hermite corrector. Integration errors and step-size control are also discussed. The end result is an accurate and efficient numerical integrator for solving ordinary differential equations and determining the sensitivity coefficients of the solution.