Suppose $f:\mathbb R\to\mathbb R^n$ satisfies $$ f'(t) = A(t)f(t), $$ where $A$ is a smooth matrix-valued function. If I know that the matrix $A(t)$ is asymptotically nilpotent, how could I prove a sub-exponential estimate for the solution $f$?
To be more explicit, suppose $A(t)^2\to0$ but $A(t)\not\to0$ as $t\to\infty$. Then one would expect slower than exponential (perhaps even linear) growth for $f$; if $A$ and $A^{-1}$ had roughly constant norm, then one would expect exponential growth. My main interest is in the case when $A(t)^2\to0$, but also $A(t)^k\to0$ for $k>2$ is interesting.
If I apply Grönwall's inequality to the function $t\mapsto|f(t)|^2$ and observe that $$ \frac{d}{dt}|f(t)|^2 = 2\langle f(t),A(t)f(t)\rangle \leq 2\|A(t)\|\cdot|f(t)|^2, $$ I get the exponential estimate $$ |f(t)| \leq |f(0)|\exp\left(\int_0^t\|A(s)\|ds\right) $$ for $t>0$. This estimate is much worse than I would expect in an asymptotically nilpotent case, but I don't know how to get a polynomial (or other sub-exponential) estimate.
Example: $n=2$ and $A(t)=\begin{pmatrix}0&1\\(1+t^2)^{-2}&0\end{pmatrix}$. Now $A(t)^2=(1+t^2)^{-2}I$ which goes to zero as $t\to\infty$. The solution to our ODE with $f(0)=(a,b)$ is $$ f(t) = \begin{pmatrix} \sqrt{1+t^2}(a+b\arctan(t)) \\ \frac1{\sqrt{1+t^2}}(at+b+bt\arctan(t)) \end{pmatrix}. $$ The solution grows essentially linearly: $|f(t)|\leq C|f(0)|(1+t)$ for any $t>0$ and some constant $C$. On the other hand, if I use Grönwall's inequality, I have the estimate $$ 2\langle f(t),A(t)f(t)\rangle = 2f_1(t)f_2(t)[1+(1+t^2)^{-2}] \leq |f(t)|^2[1+(1+t^2)^{-2}], $$ which cannot be significantly improved. Plugging this into Grönwall's inequality gives an exponential growth estimate for $f$, which much weaker than the linear estimate from the explicit solution. [The example ends here.]
I could promote the ODE to a second order one: $f''(t)=[A(t)^2+A'(t)]f(t)$. Now the coefficient $A(t)^2$ is asymptotically small, but $A'(t)$ need not be. And even if it were, I don't know how to use Grönwall for a second order ODE. If $A$ was constant, I could use nilpotency to get $f(t)=e^{At}f(0)=(I+At)f(0)$. There is a series expansion also for time-dependent $A$ (the Dyson series), but I couldn't see how to turn that into a rigorous estimate. I do not assume that $A(t)$ is nilpotent for any $t$, just that some power tends to zero as $t\to\infty$.
Question: Given some assumptions on the decay rate of $A(t)^2$ (or $A(t)^k$ for some $k>2$), what tools could I use to prove a growth estimate for norm of the solution $f(t)$? I am looking for an estimate that I could play with to see how different decay rates for $A^2$ give different growth rates for $f$.
Edit: If we denote $B(t)=A(t)+\phi(t)I$ for some scalar function $\phi$ and $g(t)=\exp\left(\int_0^t\phi(s)ds\right)f(t)$, then $g'(t)=B(t)g(t)$. One could try to get estimates for $g$ and convert them to estimates for $f$, but it seems to me that this method cannot add much. (The exponentials of integrals coming from this change of functions and Grönwall's estimate cancel each other.) This is a generalization of an idea Normal Human gave in a comment below (there $\phi$ was constant).
There is a paper by P. Hartman and A. Wintner Asymptotic Integrations of Linear Differential Equations, which proves that if $y$ satisfies the differential equation $$ y'(t) = (J + G(t))y(t), $$ where $J$ is a $n\times n$ matrix of complex numbers with eigenvalues $\lambda_1,\ldots,\lambda_n$, and $G(t)$ is a matrix of continuous complex-valued functions on $t\geq 0$ such that $|G(t)|\to0$ as $t\to\infty$, then there exist $n$ linearly independent solutions $y_1,\ldots,y_n$, satisfying the asymptotic bound $$ \log\|y_j(t)\| = (\Re\lambda_j)t + o(t), \qquad t\to\infty. $$
In your case, setting $J = \lim_{t\to\infty}A(t)$, $G(t)=A(t)-J$, the result is that $\|y(t)\| = e^{o(t)}$, because all the eigenvalues of the nilpotent matrix $J$ are zero. They also say that in the case when all $\Re\lambda_j$ are identical, it is sufficient that $$ \frac{1}{T}\int_0^T G(t)\,dt\to0. $$
They say this result was originally proved by O. Perron in Über Stabilität und asymptotisches Verhalten der Integrale von Differentialgleichungssystemen.
Regarding polynomial bounds, they say that when $$ \int_0^\infty t^{2(n-1)}\|G(t)\|\,dt < \infty $$ there is one solution $y(t)$ of the full ODE $y'=Ay$ for each solution $y_0(t)$ of the unperturbed ODE $y_0'=Jy_0$, satisfying (for a nilpotent $J$) $$ y(t) - y_0(t) = o(1), $$ which covers your example.
Perhaps some more modern textbooks on ODEs (esp. perturbation theory, structural stability) would have clearer expositions — the paper is about 40 pages long, and I didn't check the details, only the theorem statements.