Spectral unmixing is an important task in hyperspectral data exploitation. It amounts to estimating the abundance of pure spectral constituents (endmembers) in each (possibly mixed) observation collected by the imaging instrument. In recent years, several endmember extraction algorithms (EEAs) have been proposed for automated endmember extraction from hyperspectral data sets. Traditionally, EEAs extract/select only one single standard endmember spectrum for each of the presented endmember classes or scene components. The use of fixed endmember spectra, however, is a simplification since in many cases the conditions of the scene components are spatially and temporally variable. As a result, variation in endmember spectral signatures is not always accounted for and, hence, spectral unmixing can lead to poor accuracy of the estimated endmember fractions. Here, we address this issue by developing a simple strategy to adapt available EEAs to select multiple endmembers (or bundles) per scene component. We run the EEAs in randomly selected subsets of the original hyperspectral image, and group the extracted samples of pure materials in a bundle using a clustering technique. The output is a spectral library of pure materials, extracted automatically from the input scene. The proposed technique is applied to several common EEAs and combined with an endmember variability reduction technique for unmixing purposes. Experiments with both simulated and real hyperspectral data sets indicate that the proposed strategy can significantly improve fractional abundance estimations by accounting for endmember variability in the original hyperspectral data.