Methylene blue [3,7-Bis(di-methylamino) phenothiazin-5-ium chloride] is a phenothiazine dye with applications as a sensitizer for photodynamic therapy, photoantimicrobials, and dye-sensitized solar cells. Time-dependent density functional theory (TDDFT), based on (semi)local and global hybrid exchange-correlation functionals, fails to correctly describe its spectral features due to known limitations for describing optical excitations of π-conjugated systems. Here, we use TDDFT with a non-empirical optimally tuned range-separated hybrid functional to explore the optical excitations of gas phase and solvated methylene blue. We compute solvated configurations using molecular dynamics and an iterative procedure to account for explicit solute polarization. We rationalize and validate that by extrapolating the optimized range separation parameter to an infinite amount of solvating molecules, the optical gap of methylene blue is well described. Moreover, this method allows us to resolve contributions from solvent-solute intermolecular interactions and dielectric screening. We validate our results by comparing them to first-principles calculations based on the GW+Bethe-Salpeter equation approach and experiment. Vibronic calculations using TDDFT and the generating function method account for the spectra's subbands and bring the computed transition energies to within 0.15 eV of the experimental data. This methodology is expected to perform equivalently well for describing solvated spectra of π-conjugated systems.