How to estimate solutions to an ODE with an asymptotically nilpotent coefficient? 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).
 A: 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.
