We demonstrate how the functional derivatives appearing in perturbative timedependent density functional theory can be calculated using automatic differentiation. The approach starts from a computer implementation of the exchange-correlation energy functional, from which arbitrary-order derivatives are generated automatically. Automatic differentiation is shown to provide an accurate, general, and efficient implementation of higher-order exchangecorrelation functional derivatives that is easy to maintain. When used in combination with an arbitrary-order response solver, the methodology allows us to generate arbitrary-order response functions from time-dependent density functional theory. © 2010 American Chemical Society.