system of implicit nonlinear differential equations Here I have a system of nonlinear differential equations:
$
(M+2m)\ddot{x} + m(l_1 \ddot{\theta}_1\cos\theta_1 - l_1\dot{\theta}_1^2\sin\theta_1) + m(l_2\ddot{\theta}_2\cos\theta_2-l_2\dot{\theta}_2^2\sin\theta_2) = F
$
$
l_1\ddot{\theta}_1 + \ddot{x}\cos\theta_1 - g\sin\theta_1 = 0
$
$
l_2\ddot{\theta}_2 + \ddot{x}\cos\theta_2 - g\sin\theta_2 = 0
$
States are defined by me like this:
$x_1 = \theta_1 , x_2 = \dot{x_1} = \dot{\theta}_1,x_3 = \theta_2, x_4 = \dot{x_3} = \dot{\theta}_2,x_5 = x, x_6 = \dot{x_5} = \dot{x}$
I mean actually what I need to do is to linearize this system about all the states are equal to $0$. But I cannot find the $\cdots$ places below
$
\dot{x_1} = x_2 \\
\dot{x_2} = \cdots \\
\dot{x_3} = x_4 \\
\dot{x_4} = \cdots \\
\dot{x_5} = x_6 \\
\dot{x_6} = \cdots \\
$
How can I find $\dot{x_2},\dot{x_4},\dot{x_6}$?
Thanks
 A: First, put your system 
$$
(M+2m)\ddot{x} + m(l_1 \ddot{\theta}_1\cos\theta_1 - l_1\dot{\theta}_1^2\sin\theta_1) + m(l_2\ddot{\theta}_2\cos\theta_2-l_2\dot{\theta}_2^2\sin\theta_2) = F
$$
$$
l_1\ddot{\theta}_1 + \ddot{x}\cos\theta_1 - g\sin\theta_1 = 0
$$
$$
l_2\ddot{\theta}_2 + \ddot{x}\cos\theta_2 - g\sin\theta_2 = 0
$$
into a matrix formulation 
$$
\underbrace{\begin{bmatrix}
M+2m & m l_1\cos\theta_1 & ml_2\cos \theta_2\\ 
\cos \theta_1 & l_1 & 0\\
\cos \theta_2 & 0 & l_2
\end{bmatrix}}_{\boldsymbol{M}(\theta_1,\theta_2)}
\dfrac{d}{dt}
\begin{bmatrix}
\dot{x}\\
\dot{\theta_1}\\
\dot{\theta_2}\\
\end{bmatrix}=
\underbrace{\begin{bmatrix}ml_1\sin \theta_1\dot{\theta}^2_1+ml_2\sin\theta_2\dot{\theta}^2_2\\g\sin\theta_1\\g\sin\theta_2 \end{bmatrix}}_{\boldsymbol{f}(\theta_1,\theta_2,\dot{\theta_1},\dot{\theta_2})}+
\begin{bmatrix}1\\0\\0 \end{bmatrix}F.
$$
by inverting the mass matrix $\boldsymbol{M}(\theta_1,\theta_2)$, this can be done with many computer algebra systems (e.g. Maple, Python, Mathematica, MATLAB symbolic toolbox) we obtain:
$$\dfrac{d}{dt}x = \dot{x}$$
$$\dfrac{d}{dt}\theta_1 = \dot{\theta}_1$$
$$\dfrac{d}{dt}\theta_2 = \dot{\theta}_2$$
$$\dfrac{d}{dt}\begin{bmatrix}\dot{x}\\\dot{\theta}_1\\\dot{\theta}_2\end{bmatrix}=\boldsymbol{M}(\theta_1,\theta_2)^{-1}\boldsymbol{f}(\theta_1,\theta_2,\dot{\theta_1},\dot{\theta_2})+\boldsymbol{M}(\theta_1,\theta_2)^{-1}\begin{bmatrix}1\\0\\0 \end{bmatrix}F$$
