In a mathematical physical problem, I came across the following linear system of differential equations (obtained upon using Fourier series expansion): \begin{align} \frac{\mathrm{d} \rho_n}{\mathrm{d} t} + \alpha \sum_{i = 1}^{\infty} \frac{\mathrm{d} \rho_i}{\mathrm{d} t} &= - H_n \bigg( H_n \, \rho_n + \frac{\phi_n}{2} \bigg) + 1 \, , \\ \frac{\mathrm{d} \phi_n}{\mathrm{d} t} &= - H_n \, \rho_n - \phi_n \, , \end{align} where $$ H_n = 2n-1 \, , \quad \alpha \in \mathbb{R} \, , \quad \text{and} \,\, 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!