The multi-configuration time-dependent Hartree (MCTDH) method was applied to study nuclear dynamics following transitions to a manifold of vibronically coupled ground 2A1 and excited 2A2 and 2B1 states of the pyridine radical cation (PRC). These states originate from ionization out of the highest occupied orbitals of pyridine, 7a1 (nσ), 1a2 (π), and 2b1 (π), respectively, and give rise to the lowest two photoelectron bands. We focus on various theoretical and computational aspects of the MCTDH method and methodology to calculate the spectrum, taking our study of the vibronically interacting 2A1, 2A2, and 2B1 states of PRC as an example. In particular, the choice of the single-particle functions (SPFs) and schemes to combine vibrational modes are discussed.