In the lectures I've heard linear response theory was introduced multiple times, however I can't tell how this is any different from solving an inhomogeneous PDE with a greens function and there's nothing linear about that to my knowledge.
Typically it's introduced for a damped harmonic oscillator with an external force $$ \ddot{x} + \gamma \dot{x} + \omega_0^2 x = F(t) $$ The Fourier transform then gives you the susceptibility $\hat{\chi}(\omega)$ so the whole equation corresponds to multiplication in the frequency domain and we can then retrieve the solution by doing the inverse transform or doing a convolution in the time domain $$ x(t) = \int\mathrm{d}t' \chi(t-t') F(t') $$ This is how I would solve the harmonic oscillator using Fourier transforms/Greens functions (since the inverse transformed susceptibility is just the corresponding Greens function) and here I haven't made any assumptions about $F(t)$. But in my 'linear response' lectures they always mention the fact that we neglect 'higher orders of $F$' and they write something like $$ x(t) = \int\mathrm{d}t' \chi(t-t') F(t') + \mathcal{O}(F^2) $$ But the $\mathcal{O}(F^2)$ should simply always be $0$ by the derivation above, so why is there such an emphasis on the 'linear response'?
To clarify what I meant by introducing the susceptibility: take the Fourier transform of the inhomgeneous ODE \begin{align*} \hat{F} &= -\omega^2 \hat{x} -i\gamma\omega \hat{x} + \omega_0^2 \hat{x} = \hat{x}(-\omega^2 - i\gamma\omega + \omega_0^2) \\ \hat{x}(\omega) &= \hat{F}(\omega) \frac{1}{\omega_0^2 - \omega^2 - i\gamma\omega} = \hat{F}(\omega) \hat{\chi}(\omega) \end{align*} To be this procedure seems fairly general, I didn't assume anything about $F(t)$, just that it's some ODE and not a PDE? I'm not familiar with the 'Volterra Expansion', but I can't really see how more terms would get added to this expression that are non linear in $F$ somehow.