Many important engineering phenomena such as turbulent flow, fluid-structure interactions, and climate diagnostics are chaotic and sensitivity analysis of such systems is a challenging problem. Computational methods have been proposed to accurately and efficiently estimate the sensitivity analysis of these systems which is of great scientific and engineering interest. In this thesis, a new approach is applied to compute the direct and adjoint sensitivities of time-averaged quantities defined from the chaotic response of the Lorenz system and the double pendulum system. A stabilized time-integrator with adaptive time-step control is used to maintain stability of the sensitivity calculations. A study of convergence of a quantity of interest and its square is presented. Results show that the approach computes accurate sensitivity values with a computational cost that is multiple orders-of-magnitude lower than competing approaches based on least-squares-shadowing approach.