Aim:
In order to simulate a circuit model, I need to numerically solve a equation set which contains multiple second-order ODEs.
Equations of circuit:
The original equation set: $$ \begin{align} \ddot \phi_1 & = \frac{1}{C_{J1}}\cdot(\frac{\phi_4-\phi_1-\Phi_{a1}}{L_1}-I_{c1}\sin\phi_1-\frac{\dot\phi_1}{R_{J1}}) \\ \ddot \phi_2 & = \frac{1}{C_{J2}}\cdot(\frac{\phi_4-\phi_2}{L_2}-I_{c2}\sin\phi_2-\frac{\dot\phi_2}{R_{J2}}) \\ \ddot \phi_4 & = \frac{1}{C_{d1}}\cdot(\frac{\phi_5-\phi_4}{L_{d1}}-\frac{\phi_4-\phi_1-\Phi_{a1}}{L_1}-\frac{\phi_4-\phi_2}{L_2}) \\ 0 & = \frac{\phi_9-\phi_6-\Phi_{a2}}{L_3}+\frac{\phi_9-\phi_7}{L_4}-\frac{\phi_5-\phi_4}{L_{d1}} \\ \ddot\phi_6 -\ddot\phi_5 &= \frac{1}{C_{J3}}\cdot\left[\frac{\phi_9-\phi_6-\Phi_{a2}}{L_3}-I_{c3}\sin(\phi_6-\phi_5)-\frac{\dot\phi_6-\dot\phi_5}{R_{J3}}\right] \\ \ddot\phi_7 -\ddot\phi_5 &= \frac{1}{C_{J4}}\cdot\left[\frac{\phi_9-\phi_7}{L_4}-I_{c4}\sin(\phi_7-\phi_5)-\frac{\dot\phi_7-\dot\phi_5}{R_{J4}}\right] \\ 0 & =\frac{\phi_9-\phi_6-\Phi_{a2}}{L_3}+\frac{\phi_9-\phi_7}{L_4}-i_b \end{align} $$ $I_x, R_x, C_x,L_x, i_b, \Phi_{x}$ are all known constant quantities. Only phase $\phi_n$ is a time dependent variable. $\dot \phi_n$ and $\ddot \phi_n$ are first-order and second-order time derivative of phase $\phi_n$.
The simplified equation set: $$ \begin{align} \ddot\phi_1 &= f(\phi_4, \dot\phi_1, \phi_1,t)\\ \ddot\phi_2 &= f(\phi_4, \dot\phi_2, \phi_2,t)\\ \ddot\phi_4 &= f(\phi_5, \phi_4, \phi_2, \phi_1,t)\\ 0 &= f(\phi_9, \phi_7, \phi_6, \phi_5, \phi_4,t)\\ \ddot \phi_6-\ddot\phi_5 &= f(\phi_9-\phi_6,\dot\phi_6-\dot\phi_5, \phi_6-\phi_5,t)\\ \ddot \phi_7-\ddot\phi_5 &=f(\phi_9-\phi_7, \dot\phi_7-\dot\phi_5, \phi_7-\phi_5,t)\\ 0 &= f(\phi_9,\phi_7, \phi_6,t) \end{align} $$ There are total 7 equations and 7 unknowns($\phi_1, \phi_2, \phi_4, \phi_5, \phi_6, \phi_7, \phi_9$), so it must have a solution. I need to numerically solve these 7 equations simultaneously.
There are total 6 second-order items( $\ddot \phi_1,\ddot \phi_2,\ddot \phi_4,\ddot \phi_5,\ddot \phi_6,\ddot \phi_7 $ ), in order to numerically solve it, I need to set 12 initial values for them. The initial values are: $$ \dot \phi_1=\dot \phi_2=\dot \phi_4=\dot \phi_5=\dot \phi_6=\dot \phi_7=0\\ \phi_1= \phi_2= \phi_4= \phi_5= \phi_6= \phi_7=0 $$ ps: I am not sure if I set the initial conditions right or wrong, maybe I should set initial conditions for other items.
Notice that the LHS of the fifth and sixth equation is ($\ddot \phi_6-\ddot\phi_5$) and ($\ddot \phi_7-\ddot\phi_5$) , here comes the problem, I don't know how to numerically solve them. I can numerically solve the General second-order ODE, like $$ \ddot y = f(\dot y, y, t) $$ but I don't know how to solve my multi-variables coupling equation set, like ($\ddot \phi_6-\ddot\phi_5$) and ($\ddot \phi_7-\ddot\phi_5$) in it.
Could anyone teach me how to transform my equations, so I can then use some ODE solvers (like ode45 in matlab) to solve this circuit equations.
Regards,
Mengfei