Time-dependent density functional theory provides a first principles method for the calculation of frequency-dependent polarizabilities, hyperpolarizabilities, excitation energies and many related response properties. In recent years, the molecular results obtained by several groups have shown that this approach is in general more accurate than the time-dependent Hartree-Fock approach, and is often competitive in accuracy with computationally more demanding conventional ab initio approaches. In this paper, our implementation of the relevant equations in the Amsterdam Density Functional program is described. We will focus on certain aspects of the implementation which are necessary for an efficient evaluation of the desired properties, enabling the treatment of large molecules. Such an efficient implementation is obtained by: using the full symmetry of the molecule, using a set of auxiliary functions for fitting the (zeroth- and first-order) electron density, using a highly vectorized and parallelized code, using linear scaling techniques, and, most importantly, by solving the response equations iteratively.