First-principle simulations aimed at accurately reproducing the excited state properties of a large series of ladder-type π-conjugated organic molecules containing heteroatoms (Si, S, B, O, and N) have been performed. In particular, time-dependent density functional theory (TD-DFT) calculations relying on several global and range-separated hybrid functionals have been carried out in conjunction with three variations of the polarizable continuum model (PCM), namely, the linear-response (LR), corrected linear-response (cLR), and state-specific (SS) approaches. For this series of molecules, similar to many borate derivatives, the cLR-PCM-TD-M06-2X approach can be used to reproduce the auxochromic effects that tune the 0-0 energies. However, TD-DFT yields rather large absolute deviations with respect to the experimental 0-0 energies. These systematic errors can be reduced by more than 0.1 eV when scaled opposite spin-configuration interaction singles with a double correction [SOS-CIS(D)] vertical calculations are combined to the PCM-TD-DFT results. This study demonstrates that such a "hybrid" scheme, where the geometrical and vibrational parameters, as well as the solvation effects, are determined with PCM-TD-DFT, whereas the transition energies are obtained with a wavefunction-based method, offers a useful compromise between accuracy and computational cost.