An accurate determination of frequency-dependent molecular hyperpolarizabilities is at the same time of possible technological importance and theoretically challenging. For large molecules, Hartree-Fock theory was until recently the only available ab initio approach. However, correlation effects are usually very important for this property, which makes it desirable to have a computationally efficient approach in which those effects are (approximately) taken into account. We have recently shown that frequency-dependent hyperpolarizabilities can be efficiently obtained using time-dependent density functional theory. Here, we shall present the necessary theoretical framework and the details of our implementation in the Amsterdam Density Functional program. Special attention will be paid to the use of fit functions for the density and to numerical integration, which are typical of density functional codes. Numerical examples for He, CO, and para-nitroaniline are presented, as evidence for the correctness of the equations and the implementation.