We present a discrete solvent reaction field (DRF) model for the calculation of frequency-dependent hyperpolarizabilities of molecules in solution. In this model the solute is described using density functional theory (DFT) and the discrete solvent molecules are described with a classical polarizable model. The first hyperpolarizability is obtained in an efficient way using time-dependent DFT and the (2n+1) rule. The method was tested for liquid water using a model in which a water molecule is embedded in a cluster of 127 classical water molecules. The frequency-dependent first and second hyperpolarizabilities related to the electric field induced second harmonic generation (EFISH) experiment, were calculated both in the gas phase and in the liquid phase. For water in the gas phase, results are obtained in good agreement with correlated wave function methods and experiments by using the so-called shape-corrected exchange correlation (xc)-potentials. In the liquid phase the effect of using asymptotically correct functionals is discussed. The model reproduced the experimentally observed sign change in the first hyperpolarizaibility when going from the gas phase to the liquid phase. Furthermore, it is shown that the first hyperpolarizability is more sensitive to damping of the solvent-solute interactions at short range than the second hyperpolarizability. (C) 2003 American Institute of Physics.