I have this system of differential equations which describes the motion of a missile launcher model with 5 degrees of freedom:
(1)$$(m_w +m_v +m_p)\ddot{y}_w - (m_v + m_p)h_v\ddot{\vartheta}_w \sin(\vartheta_w) + m_p \xi_p \ddot{\vartheta}_v \cos(\vartheta_v) \\ - (m_v + m_p)h_v\dot{\vartheta}^2_w \cos(\vartheta_w) - m_p\xi_p \dot{\vartheta}^2_v \sin(\vartheta_v) + m_p \dot{\xi}_p \dot{\vartheta}_v \cos(\vartheta_v) +k_{w11}\lambda_{w11} \\ + k_{w12}\lambda_{w12} -c_{w11}\dot{\lambda}_{w11} -c_{w12}\dot{\lambda}_{w12} = -(m_w +m_v +m_p)g$$
(2)$$[I_w+I_{p_{zp}}+I_v+(m_v+m_p)h^2_v]\ddot{\vartheta}_w + [I_{p_{zp}}+m_p h_v \xi_p(\cos\vartheta_w\sin\vartheta_v \\ - \sin\vartheta_w\cos\vartheta_v)]\ddot\vartheta_v - (m_v + m_p)h_v\ddot{y}_w \sin\vartheta_w - m_p h_v \ddot{\xi}_p(\cos\vartheta_w\cos\vartheta_v \\ + \sin\vartheta_w\sin\vartheta_v) + 2m_p h_v \dot{\xi}_p \dot{\vartheta}_w \sin\vartheta_w\cos\vartheta_v \\ -2m_p h_v \dot{\xi}_p \dot{\vartheta}_v \sin\vartheta_w\cos\vartheta_v +m_p h_v \xi_p \dot{\vartheta}^2_v(\cos\vartheta_w\cos\vartheta_v + \sin\vartheta_w\sin\vartheta_v)\\ +k_{w11}l_{w1}\lambda_{w11} + k_{w12}l_{w2}\lambda_{w12} - k_{w31}\lambda_{w31} -c_{w11}l_{w1}\dot{\lambda}_{w11}+c_{w12}l_{w2}\dot{\lambda}_{w12} \\ + c_{w31}\dot{\lambda}_{w31} = (m_v + m_p)gh_v\sin(\vartheta_w+\vartheta_{wst})$$
(3)$$(I_v + I_{p_{zp}} + m_p\xi^2_p)\ddot{\vartheta}_v + m_p \xi_p \ddot{y}_w\cos\vartheta_v + [I_v + I_{p_{zp}} + m_p h_v \xi_p(\cos\vartheta_w\sin\vartheta_v \\ - \sin\vartheta_w\cos\vartheta_v)]\ddot{\vartheta}_w + 2m_p\xi_p\dot{\xi}_p\dot{\vartheta}_v + m_p h_v\xi_p\dot{\vartheta}^2_w(\sin\vartheta_w\sin\vartheta_v \\ + \cos\vartheta_w\cos\vartheta_v)+k_{w31}\lambda_{w31}-c_{w31}\dot{\lambda}_{w31} = -m_pg\xi_p\cos(\vartheta_v + \vartheta_{vst})$$
(4)$$ m_p\ddot{\xi}_p-m_ph_v\ddot{\vartheta}_v(\cos\vartheta_w\cos\vartheta_v + \sin\vartheta_w\sin\vartheta_v) \\ + m_p\ddot{y}_w\sin\vartheta_w + m_ph_v\dot{\vartheta}^2_w(\sin\vartheta_w\cos\vartheta_v - \cos\vartheta_w\sin\vartheta_v) \\ + m_p\xi_p\dot{\vartheta}^2_v = -m_pg\sin(\vartheta_v + \vartheta_{vst}) + P_p$$
(5)$$I_{p_{xp}}\ddot{\varphi}_p = M_p$$.
The scheme of the model is this: scheme of the launcher.
I know all the coefficients but I want to see the trend in the first 2 seconds of the variables $y_w, \vartheta_w, \vartheta_v, \xi_p, \phi_p$. Unfortunately I tried with Matlab ode45 but the system is coupled and nonlinear, can you suggest me a numerical method, a code, or anything that may help me solve this system?
After understanding the procedure I can do some more advanced calculations but I need to figure out the approach to solve it first. Thanks, really thanks, to anybody who wants to help.
Loris
EDIT: same equations but with numbers and with linearised sines and cosines so that you can clearly see where is the difficulty:
(1.1)$$\ddot{y}_w = 0.0917\ddot{\vartheta_w}\vartheta_w -0.1667\xi_p\ddot{\vartheta_w} +0.0917\dot{\vartheta_w}^2 + 0.1667\xi_p\dot{\vartheta_v}^2\vartheta_v -0.1667\dot{\xi_p}\dot{\vartheta_v} -4166.67(2y_w +0.3916\vartheta_w)+ 4.1667*(2\dot{y}_w +0.3916\dot{\vartheta_w}) + 9.81$$
(2.1)$$\ddot{\vartheta_w} = -0.3217\ddot{\vartheta_v} -0.4531\xi_p\ddot{\vartheta_v}(\vartheta_v-\vartheta_w) +1.2082\ddot{y}_w\vartheta_w +0.4531\ddot{\xi}_p(1+\vartheta_w\vartheta_v) -0.9602\dot{\xi}_p\vartheta_w(\vartheta_w-\vartheta_v) -0.4531\xi_p\dot{\vartheta_v}^2(1+\vartheta_w\vartheta_v) -4.393\cdot 10^4y_w -1.72\cdot 10^4\vartheta_w +366(\vartheta_v-\vartheta_w) +21.5\dot{y}_w +21.78\dot{\vartheta_w} +36.6(\dot{\vartheta_w}-\dot{\vartheta_v}) +11.8529\vartheta_w$$
(3.1)$$\ddot{\vartheta_v} = -3.9734(\xi_p^2\ddot{\vartheta_v}+\xi_p\ddot{y}_w) -\ddot{\vartheta_w} -0.8197\xi_p(\vartheta_v-\vartheta_w)\ddot{\vartheta_w} -7.9468\xi_p\dot{\xi}_p\dot{\vartheta_v} -0.8197\xi_p\dot{\vartheta_w}^2(1+\vartheta_w\vartheta_v) -662.23(\vartheta_v-\vartheta_w) -66.223(\dot{\vartheta_v}-\dot{\vartheta_w}) +38.98\xi_p$$
(4.1) $$\ddot{\xi}_p = 0.2063\ddot{\vartheta_v}(1+\vartheta_w\vartheta_v) -\ddot{y}_w\vartheta_w -0.2063\dot{\vartheta_w}^2(\vartheta_w-\vartheta_v) -\xi_p\dot{\vartheta_v}^2 -9.81\vartheta_v +333$$
(5.1) $$\ddot{\phi}_p = 81$$
Letting the four unknowns form the vector $\mathbf z$, your system can be written
$$\ddot{\mathbf z}=P(\dot{\mathbf z}, \mathbf z)\ddot{\mathbf z}+\mathbf q(\dot{\mathbf z}, \mathbf z)$$ where $P$ is a matrix and $\mathbf q$ a vector, both depending on $\dot{\mathbf z}$ and $\mathbf z$.
Solve this by
$$\ddot{\mathbf z}=(I-P(\dot{\mathbf z}, \mathbf z))^{-1}\mathbf q(\dot{\mathbf z}, \mathbf z).$$
Then with $\mathbb w=\dot{\mathbb z}$, you get a nonlinear system of the first order,
$$\begin{cases} \dot{\mathbf w}=(I-P({\mathbf w}, \mathbf z))^{-1}\mathbf q({\mathbf w}, \mathbf z)\\\dot{\mathbf z}=\mathbf w \end{cases}$$ ready for numerical resolution.