The response of the density matrix to an external field is calculated in the adiabatic time-dependent density functional (TDDFT) theory by mapping the equation of motion for the driven single-electron density matrix into the dynamics of coupled harmonic oscillators. The resulting nonlinear response functions and the closed expressions for nonlinear frequency-dependent polarizabilities are derived. These expressions include transition densities and frequencies calculated in the linear response TDDFT, and higher order functional derivatives of the exchange-correlation functional. Limitations of the applicability of the traditional sum over states approach for computing the nonlinear response to the TDDFT are discussed. (C) 2003 American Institute of Physics.