The "hierarchical equations of motion"(HEOM) method is a powerful exact numerical approach to solve the dynamics and find the steady-state of a quantum system coupled to a non-Markovian and nonperturbative environment. Originally developed in the context of physical chemistry, it has also been extended and applied to problems in solid-state physics, optics, single-molecule electronics, and biological physics. Here we present a numerical library in Python, integrated with the powerful QuTiP platform, which implements the HEOM for both bosonic and fermionic environments. We demonstrate its utility with a series of examples consisting of benchmarks against important known results and examples demonstrating insights gained with this library for this article. For the bosonic case, our results include demonstrations of how to fit arbitrary spectral densities with different approaches, and a study of the dynamics of energy transfer in the Fenna-Matthews-Olson photosynthetic complex. For the latter, we both clarify how a suitable non-Markovian environment can protect against pure dephasing, and model recent experimental results demonstrating the suppression of electronic coherence. Importantly, we show that by combining the HEOM method with the reaction coordinate method we can observe nontrivial system-environment entanglement on timescales substantially longer than electronic coherence alone. We also demonstrate results showing how the HEOM can be used to benchmark different strategies for dynamical decoupling of a system from its environment, and show that the Uhrig pulse-spacing scheme is less optimal than equally spaced pulses when the environment's spectral density is very broad. For the fermionic case, we present an integrable single-impurity example, used as a benchmark of the code, and a more complex example of an impurity strongly coupled to a single vibronic mode, with applications to single-molecule electronics.