How to solve multiple coupled differential equations by MATLAB?

Question

The differential equations are as follows

\begin{align} \dot x_i &= \gamma(\mu-r^2_i)x_i -\omega_iy_i +\varepsilon F(t) +\tau\sin(R_i-\phi_i) \tag{13} \\ \dot y_i &= \gamma(\mu-r^2_i)y_i +\omega_ix_i \tag{14} \\ \dot \omega_i &= -\varepsilon F(t)\frac{y_i}{r_i} \tag{15} \\ \dot \alpha_i &= \eta x_i F(t) \tag{16} \\ \dot \phi_0 &= 0 \tag{17} \\ \dot \phi_i &= \sin\left(R_i-\operatorname{sgn}(x_i)\cos^{-1}\left(-\frac{y_i}{r_i}\right)-\phi_i\right), \quad\forall\ i\neq 0 \tag{18} \\ \end{align} with $$R_i = \frac{\omega_i}{\omega_0} \operatorname{sgn}(x_0) \cos^{-1}\left(-\frac{y_0}{\sqrt{x^2_0+y^2_0}}\right) \tag{19}$$ and $$F(t) = P_{\text{teach}}(t) -\sum_{i=0}^N \alpha_i x_i \tag{20}$$ $$P_{\text{teach}}=0.8 \sin(15t) +\cos(30t)-1.4 \sin(45t)-0.5 \cos(60t) \tag{21}$$ And the initial conditions are as follows $$\alpha_i(0) = \phi_i(0) = 0, x_i(0) = 1, y_i(0) = 0\ \forall\ i, \mu=1, \gamma=8, \varepsilon=0.9, \eta=0.5, \tau=2.$$ The frequencies $\omega_i(0)$ are uniformly distributed from $6$ to $70$.

---

Based on the image description, we have $5$ variables, $x$, $y$, $\omega$, $\alpha$ and $\phi$.

So how to solve the ODEs if we have $4$ equation sets ($i = 1, 2, 3, 4$)? And the $x$, $y$, $\omega$, $\alpha$ and $\phi$ will be vectors respectively.

I just thought that the ode45 command can solve it. In that case, we must split the ODEs to $4 \times 5=20$ equations?

However, if we get many equation sets (for example $i= 1, 2, 3, \ldots, 100$), how should we solve this problem?

Differential equations are derived from this paper:

Righetti L, Buchli J, Ijspeert A J.: From Dynamic Hebbian Learning for Oscillators to Adaptive Central Pattern Generators[C] || Proceedings of $3$rd International Symposium on Adaptive Motion in Animals and Machines -- AMAM $2005$. Verlag ISLE, Ilmenau, $2005$.

Answer 1

You should use the vector and elementwise array operations that Matlab provides. That is, you take the input vector and cut it into the parts for $x,y,...$ and then do the operations on them never even declaring an index. Something like

dotu = odefunc(t,u)
    x = u(1:N)
    y = u(N+1,2*N)
    ....

    r2 = x.^2+y.^2
    r = r2.^0.5
    ....
    doty = gamma*(mu-r2).*y+omega.*x
    Ft=F(t)
    dotomega = -epsilon*Ft*y./r
    ...

and in the end join the derivatives back together in the order of the state vector.