After thinking out how to convert a non-homogeneous linear differential equation, with a polynomial input, to a homogeneous linear differential equation in general for this question I started playing around with this myself for fun.
I used a driven damped harmonic oscillator and used a Taylor polynomial, $p_{in}(t)$, of a sine as input
$$ \begin{bmatrix} \dot{x} \\ \ddot{x} \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ -k & -d \end{bmatrix} \begin{bmatrix} x \\ \dot{x} \end{bmatrix} + p_{in}(t) \begin{bmatrix} 0 \\ 1 \end{bmatrix}. $$
By using the method described in my answer to the linked question (not sure if this already has a name in the literature) this can be transformed in to a homogeneous equation in $y(t)$,
$$ x(t) = y(t) + p_{out}(t), $$
such the polynomial, $p_{out}(t)$, is constructed such that,
$$ \begin{bmatrix} \dot{y} \\ \ddot{y} \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ -k & -d \end{bmatrix} \begin{bmatrix} y \\ \dot{y} \end{bmatrix}. $$
My hypothesis was that $y(t)$ would act as the transient response and $p_{out}(t)$ as the steady state response of $x(t)$. For a sine wave of magnitude one and frequency $\omega$ the steady state response of this system, using Fourier transform, can be found to be,
$$ x_{ss}(t) = \frac{k-\omega^2}{(k-\omega^2)^2+\omega^2d^2}\sin{\omega t} - \frac{\omega d}{(k-\omega^2)^2+\omega^2d^2}\cos{\omega t}. $$
Instead of using an actual sine as input, but the Taylor polynomial at $t=0$ I suspected that $p_{out}(t)$ would approximate $x_{ss}(t)$ near $t=0$ as well. This appeared to be correct under certain circumstances, namely when $\omega$ is below the resonant frequency of the system, defined as $\sqrt{k-\frac{d^2}{4}}$. When $\omega$ is above the resonant frequency, then $p_{out}(t)$ seems to be dominated by the transient response near $t=0$, one example can be seen in the following figure,
where the the blue line is equal to $x_{ss}(t)$ and the red dotted line is equal to $p_{out}(t)$, which is found from $p_{in}(t)$, which is a Taylor polynomial approximation of order 91. For the used system I used $k=1$, $d=0.8$ and the forcing frequency used to find the coefficients of $p_{in}(t)$ is $1.031$ rad/s or $1.125$ times the resonant frequency.
Below the resonant frequency there will also be some transient response present, but does not dominate over the steady state response. So what is so special about the resonant frequency in this case?
Edit:
After thinking about this a bit more the reason why any transient response is present can be assigned to the fact that any polynomial will be approximately be equal to its highest power when $|t|>>0$ and thus will go of to $\pm\infty$. So when the error between the input polynomial and the sine wave drops significantly, the output still has to catch up with an exponential decay from the offset from zero. Lower the forcing frequency, $\omega$, will move the moment, at which the error drops significantly, back in time, since the coefficient of the highest power, $m$, of a Taylor polynomial of a sine wave is proportional to $\omega^m$. But the exponential decay rate will be constant for the same system/matrix, thus lowering the forcing frequency will lower the amplitude of a transient response compared to the steady state response near $t=0$.
Before the resonant frequency the amplitude of the steady state response will remain roughly constant, near it the amplitude might increase (depending on whether the system under-, over- or critically damped), however above it the amplitude will drop of with $\propto\omega^{-2}$. Thus when the forcing frequency is above the resonant frequency the amplitude of the steady state response will become smaller relative to the transient response.
Both these processes should increase the amplitude of a transient like response relative to the steady state response, but I do not know how to assign actual numerical values to these processes for a given system and forcing frequency.