Algorithms are presented for overcoming the computational challenge of nonlinear response functions which describe the response of a classical system to a sequence of n pulses and depend on nth order multipoint stability matrices containing signatures of chaos. Simulations for the Lorentz gas demonstrate that finite field algorithms can be effectively used for the robust, long time calculation of nonlinear response functions. These offer the possibility to characterize chaos beyond the commonly used Lyapunov exponents and suggest new experimentally accessible measures of chaos.