Time dependent self-consistent equations

I am facing the following problem. I need to solve numerically a set of coupled equations

$$i\frac{d}{dt}f_{n}^{(i)}(t) = \left[U\cdot n(n-1) + \mu\cdot n\right]f_{n}^{(i)}(t) - \sqrt{n+1}\Phi_i^{*}\ f_{n+1}^{(i)}-\sqrt{n}\Phi_{i}\ f_{n-1}^{(i)}$$ where $U,\mu$ are just constants and $$\Phi_{i} = \sum\limits_{n=1}^{N}\left( f_{n-1}^{(i+1)}(t)\right)^{*}\ f_{n}^{(i+1)}(t)\sqrt{n}$$ has to be determined self-consistently. These are coupled differential equations where $i = 1,2,\ldots, M$ and $n = 0,1,,\ldots, N$. What kind of numerical methods can be used in order to solve this problem efficiently? To be clear, I am looking for time evolution of each $f_{n}^{(i)}$ (complex numbers).

• What about Runge-Kutta algorithm for $(N+1)M$ variables? – WoofDoggy Sep 30 '15 at 14:07
• I don't see anything unusual here -- the usual methods (e.g. Runge-Kutta, linear multistep) should work fine. – David Ketcheson Sep 30 '15 at 17:17
• Have you looked at the literature on this topic? – Jeff Dec 23 '15 at 4:09
• Thanks for the comments. Indeed, there is nothing unusual in those equations and Runge-Kutta can handle it. For slowly oscillating solutions exponent method also works pretty well. – WoofDoggy Dec 23 '15 at 8:32