Vtyjkyuk

In a mathematical physical problem, I came across the following linear system of differential equations (obtained upon using Fourier series expansion):
beginalign
fracmathrmd rho_nmathrmd t +
alpha sum_i = 1^infty fracmathrmd rho_imathrmd t
&= - H_n bigg( H_n , rho_n
+ fracphi_n2 bigg) + 1 , , \
fracmathrmd phi_nmathrmd t
&= - H_n , rho_n - phi_n , ,
endalign
where
$$
H_n = 2n-1 , ,
quad
alpha in mathbbR , ,
quad
textand
,,
n ge 1, .
$$

The system is subject to the initial conditions $rho_n(0)=phi_n(0)=0$.

The goal is to determine the general expression of the coefficients $rho_n$ and $phi_n$.
When $alpha=0$, the solution of the problem is easy and straightforward.

I have tried to solve the above linear system of differential equations for $alpha ne 0$ using the Laplace transform technique but without success. The calculation of the inverse Laplace transform doesn't seem to be possible (since the Laplace-transformed function has an infinite number of singularities and choosing an abscissa of convergence $sigma>0$ for which the contour is located to the right of all singularities is not within reach.)

I was wondering whether someone here could be of help and try to tell how one can solve such a mathematical problem.
Your hints and ideas and very welcome.
Very much thanks!