A resonance phenomenon appears in the Raman response when the exciting light has frequency close to electronic transitions. Unlike for molecules and for graphene, the theoretical prediction of such frequency-dependent Raman response of crystalline systems has remained a challenge. Indeed, many Raman intensity first-principle calculations are nowadays done at vanishing light frequency, using static Density-Functional Perturbation Theory, thus neglecting the frequency dependence and excitonic effects. In this work, we propose a finite-difference method to compute frequency-dependent Raman intensities. Recently, we used this methodology for the computation of the first-order frequency-dependent Raman intensity [1], with excitonic effects described by the Bethe-Salpeter equation. We found these to be crucial for the accurate description of the experimental enhancement for laser photon energies around the gap. We have also generalized this approach to the more complex second-order Raman intensity, with phonon losses coming from the entire Brillouin zone. Interestingly, even without excitonic effects, one is able to capture the main relative changes in the frequency-dependent Raman spectrum at fixed laser frequencies. However, excitonic effects might affect significantly the intensity of specific modes and also lead to a global increase of absolute intensities. [1] Y. Gillet, M. Giantomassi, X. Gonze, Phys. Rev. B 88, 094305 (2013