Boron-dipyrromethene (BODIPY) molecules form a class of fluorescent dyes known for their exceptional photoluminescence properties. Today, they are used extensively in various applications from fluorescent imaging to optoelectronics. The ease of altering the BODIPY core has allowed scientists to synthesize dozens of analogues by exploring chemical substitutions of various kinds or by increasing the length of conjugated groups. However, predicting the impact of any chemical change accurately is still a challenge, especially as most computational methods fail to describe correctly the photophysical properties of BODIPY derivatives. In this study, the recently developed coupled cluster method called "domain-based local pair natural orbital similarity transformed equation of motion-coupled cluster singles and doubles" (DLPNO-STEOM-CCSD) is employed to compute the lowest vertical excitation energies of more than 50 BODIPY molecules. The method performs remarkably well yielding an accuracy of about 0.06 eV compared to the experimental absorption maxima. We also provide an estimate to the error made by neglecting vibronic effects in the computed spectra. The dyes selected for investigation here span a large range of molecular sizes and chemical functionalities and are embedded in solvents with different polarities. We have also investigated if the method is able to correctly reproduce the impact of a single chemical modification on the absorption energy. To characterize the method in more specific terms, we have studied four large BODIPY analogues used in real-life applications due to their interesting chemical properties. These examples should illustrate the capacity of the DLPNO-STEOM-CCSD procedure to become a method of choice for the study of photophysical properties of medium to large organic compounds.