Series expansion of ODE solution Consider the ODE for $y(t) \in \mathbb R^n$ with square matrix $A$ and bilinear operator $H: \mathbb R^n \times \mathbb R^n \rightarrow \mathbb R^n$, I am looking for a closed form series for the solution to:
$$y'(t) = A y(t) -H(y(t),y(t))$$
with initial condition $y(0)=y_0$
Clearly, 
$$y(t) = e^{At} y_0 -\int_0^t e^{A(t-s)} H(y(s),y(s)) \ ds.$$
Clearly, one can now substitute this expression again into the last term which leaves us with
$$y(t) = e^{At} y_0 -\int_0^t e^{As} H(y_0,y(s)) \ ds-\int_0^t \int_0^s e^{A(t-s)} e^{A(s-s_1)} H(H(y(s_1),y(s_1)),y(s)) \ ds_1 \ ds .$$
In principle, one should be able to find a closed form infinite series expression for this solution, i.e. write $$y(t) =\sum_{n=0}^{\infty} a_n(y_0)$$
with $\vert a_n(y_0) \vert \lesssim \vert y_0 \vert^n$ but I cannot really find a way to write it down.
Does anybody know?
EDIT: Please see the comment by LutzL below for an explanation for what is needed.
 A: The power series coefficients for $y(t)=\sum_ky_kt^k$ follow the recursion
$$
(k+1)y_{k+1}=Ay_k-\sum_{j=0}^kH(y_j,y_{k-j}).\tag1
$$
Their size is bounded by the recursive inequality
$$
(k+1)\|y_{k+1}\| \le \|A\|\,\|y_k\|+\sum_{j=0}^k\|H\|\,\|y_j\|\,\|y_{k-j}\|
\tag2
$$
in the induced operator and tensor norms. Now compare this to the scalar equation
$$
u'=au+hu^2,~~a=\|A\|,~h=\|H\|,~u(t)=\sum_ku_kt^k\tag3
$$
which has an exact solution 
$$
au(t)^{-1}+h=(au_0^{-1}+h)e^{-at}\implies u(t)=\frac{u_0e^{at}}{a(a-hu_0(e^{at}-1))}\tag4
$$
The power series coefficients of $u(t)=\sum_ku_kt^k$ follow the recursion
$$
(k+1)u_{k+1}=au_k+h\sum_{j=0}^ku_ju_{k-j}.\tag5
$$
Comparing this recursion (5) with the recursive norm inequality (2) resulting from the differential equation one concludes that if $\|y_j\|\le u_j$ for $j=0,..,k$, then also $\|y_{k+1}\|\le u_{k+1}$, that is, this holds for all $k$ if $u_0=\|y_0\|$. 
The exact solution for $u$ has a pole at $t=\rho=\frac1a\ln(1+\frac{a}{hu_0})\le\frac1{hu_0}$. Additionally one sees that $u_k$ is a polynomial in $u_0$ of degree $k+1$.
From general theory about the radius of convergence one concludes that for any $r<\rho$ there is a constant $C_r$ so that $\|y_k\|\le u_k\le C_rr^{-k-1}$. For $r=\frac12\rho$ for example one thus gets $$\|y_k\|\le C(2\|H\|\,\|y_0\|)^{k+1}.$$
