We calculated the one- (OPA) and two-photon absorption (TPA) spectra of two large water-soluble stilbene derivatives presenting TPA cross sections of about 400 GM. However, the reported experimental TPA spectra present a spectral gap region, and a theoretical study of these promising molecules seems now timely and relevant. These molecules are composed of 200 or more atoms, becoming a challenge to obtain the TPA spectra even using density functional theory at the time-dependent quadratic response formalism. Thus, both OPA and TPA were also calculated using the INDO-S semi-empirical method. We used explicit solvent molecules using the sequential-quantum mechanics/molecular mechanics to include the solvent effects. Our results show that different transitions are participating in the OPA and TPA processes and that exchange-correlation functionals, including larger Hartree-Fock contributions, provide a better description of the OPA spectra; however, the opposite trend is observed on the TPA spectra. Alternatively, INDO-S/CISD, including contributions from single and double excitations, systematically describes both OPA and TPA bands with similar shifts and better reproduces the relative intensities of the two TPA bands compared to the experimental ones. The OPA spectra are characterized by a Highest Occupied Molecular Orbital-Lowest Unoccupied Molecular Orbital (HOMO-LUMO) excitation, while the low-energy TPA band is ascribed to a single transition encompassing the (HOMO-1)-LUMO and HOMO-(LUMO+1) excitations and the high-energy one is a combination of several transitions. Thus, although more studies are required to better assess the capability of the INDO-S/CISD method in describing the TPA spectra of large molecules, our results corroborate that it is a promising alternative.