We present a method for determining quadratic response molecular properties of heterogeneously solvated molecules. The method enables the calculation of frequency-dependent hyperpolarizabilities and two-photon transition matrix elements of molecules that are (i) located at the surface of a metal or (ii) solvated alongside the surface of a metal. The metal is represented as a perfect conductor and the solvent is described by a dielectric medium. The method is able to describe the electronic structure of the molecular systems both at the uncorrelated and correlated electronic structure levels. The latter is given by the multiconfigurational self-consistent field method. We have completed sample calculations of the frequency-dependent first-order hyperpolarizability tensor and two-photon matrix elements for heterogeneously solvated CO. We find that the molecular properties depend strongly on the heterogeneous solvent configuration.