We present a method of modeling vibrational resonance Raman scattering (RRS) spectra of isolated and solvated systems with the inclusion of Franck-Condon (FC) and Herzberg-Teller (HT) effects and a full account for possible differences between the harmonic potential energy surfaces of the initial and resonant electronic states. It describes fundamentals, overtones, and combination bands and computes the RRS spectrum as a two-dimensional function of the incident and scattered frequencies. The theoretical foundations of the method are described and the differences with other currently available methodologies are outlined. Applications to the phenoxyl radical in the gas phase and indolinedimethine-malononitrile (IDMN) in acetonitrile and cyclohexane solution are reported, as well as comparisons with available experimental data.