Functional MRI (fMRI) studies typically analyze data by applying a single function – across the entire brain – to relate what is measured (blood oxygenation fluctuations) to the underlying neural activity. However, this hemodynamic response function (HRF), is known to vary considerably across brain regions in healthy individuals, and even more prominently in clinical populations (e.g., AIDS, Alzheimer’s). An improved characterization of HRF variability would improve cognitive science experimentation, effective connectivity analysis, and may be crucial for early detection of certain diseases. Here, a method is suggested for altering stimulus presentation timing during task related fMRI experiments that aims to maximize characterization of HRF variability while minimizing the number of trials required to accomplish this. To do so, d-optimality constraints are applied for sparse sampling of the HRF in the temporal domain. We first demonstrate this approach using simulated data over a range of background noise fluctuations. Using simulated data, we were able to recover HRF signal estimates with <10% sum of squared error (SSE) using 73% and 47% less stimulus events using D-optimal sampling compared to fixed or random designs respectively. We then utilized this method for designing the stimulus timing in an event-related fMRI experiment. Empirically, we were able to detect the initial dip in 53% of subjects, a part of the HRF signal that is thought to reflect oxygen usage and often obscured when using conventional experimental design paradigms.