  We propose a self-validating scheme to calculate the unbiased responses of
quantum many-body systems to external fields of arbibraty strength at any
temperature. By switching on a specified field to a thermal pure quantum state
of an isolated system, and tracking its time evolution, one can observe an
intrinsic thermalization process driven solely by many-body effects. The
transient behavior before thermalization contains rich information on excited
states, giving the linear and nonlinear response functions at all frequencies.
We uncover the necessary conditions to clarify the applicability of this
formalism, supported by a proper definition of the nonlinear response function.
The accuracy of the protocol is guaranteed by a rigorous upper bound of error
exponentially decreasing with system size, and is well implemented in the
simple ferromagnetic Heisenberg chain, whose response at high fields exhibits a
nonlinear band deformation. We further extract the characteristic features of
excitation of the spin-1/2 kagome antiferromagnet; the wavenumber-insensitive
linear responses from the possible spin liquid ground state, and the
significantly broad nonlinear peaks which should be generated from numerous
collisions of quasi-particles, that are beyond the perturbative description.
