We propose a novel strategy for global sensitivity analysis of ordinary differential equations. It is based on an error-controlled solution of the partial differential equation (PDE) that describes the evolution of the probability density function associated with the input uncertainty/variability. The density yields a more accurate estimate of the output uncertainty/variability, where not only some observables (such as mean and variance) but also structural properties (e.g., skewness, heavy tails, bi-modality) can be resolved up to a selected accuracy. For the adaptive solution of the PDE Cauchy problem we use the Rothe method with multiplicative error correction, which was originally developed for the solution of parabolic PDEs. We show that, unlike in parabolic problems, conservation properties necessitate a coupling of temporal and spatial accuracy to avoid accumulation of spatial approximation errors over time. We provide convergence conditions for the numerical scheme and suggest an implementation using approximate approximations for spatial discretization to efficiently resolve the coupling of temporal and spatial accuracy. The performance of the method is studied by means of low-dimensional case studies. The favorable properties of the spatial discretization technique suggest that this may be the starting point for an error-controlled sensitivity analysis in higher dimensions.