In this paper we develop a new method for numerically approximating sensitivities in parameter-dependent ordinary differential equations (ODEs). Our approach, intended for situations where the standard forward and adjoint sensitivity analyses become too computationally costly for practical purposes, is based on the Peano-Baker series from control theory. Using this series, we construct a representation of the sensitivity matrix S and, from this representation, a numerical method for approximating S. We prove that, under standard regularity assumptions, the error of our method scales as O(Δtmax2), where Δtmax is the largest time step used when numerically solving the ODE. We illustrate the performance of the method in several numerical experiments, taken from both the systems biology setting and more classical dynamical systems. The experiments show the sought-after improvement in running time of our method compared to the forward sensitivity approach. In experiments involving a random linear system, the forward approach requires roughly n longer computational time, where n is the dimension of the parameter space, than our proposed method.