Suppose i have a system of coupled oscillators: $\ddot{\phi}_i = \sum_j M_{ij}\phi_j $. This is equivalent to the Matrix differential equation: $$ \begin{pmatrix} \dot{\phi}_1\\ \vdots \\ \dot{\phi}_N \\ \dot{v}_1 \\ \vdots\\ \dot{v}_N\\ \end{pmatrix} = \begin{pmatrix} 0 && \mathbb{I}\\ M && 0\\ \end{pmatrix} \begin{pmatrix} \phi_1 \\ \vdots \\ \phi_N \\ v_1 \\ \vdots\\ v_N\\ \end{pmatrix} $$ where $\dot{v}_i = \ddot{\phi}_i$, right?
I know from diagonalization the eigenvalues and eigenvectors of $M$. How are they linked to the eigenvalues and eigenvectors of $\begin{pmatrix} 0 && \mathbb{I}\\ M && 0\\ \end{pmatrix}$? Therefore, how do I find the solutions $\phi_i$ based on the eigenvalues and eigenvectors of $M$?
Since you can diagonalize the symmetric matrix M, to $$M_{ij} \mapsto D_{ij}= \lambda_i \delta_{ij} \leadsto \\ D=\begin{pmatrix} \lambda_1&...&0 \\ 0& \ddots&0 \\ 0& ...&\lambda_N \end{pmatrix} ~~~~ \leadsto ~~~~ \ddot{\varphi}_i = \lambda _i \varphi_i, $$ you already have the solutions! (Negative semidefinite eigenvalues.)
In phase space, you have $\dot \varphi_i = u_i$, so that $\dot u_i=\ddot \varphi_i$ $\implies \ddot u_i= \lambda_i u_i$ , and $$ \begin{pmatrix} \dot{\varphi}_1\\ \vdots \\ \dot{\varphi}_N \\ \dot{u}_1 \\ \vdots\\ \dot{u}_N\\ \end{pmatrix} = \begin{pmatrix} 0 && \mathbb{I}\\ D && 0\\ \end{pmatrix} \begin{pmatrix} \varphi_1 \\ \vdots \\ \varphi_N \\ u_1 \\ \vdots\\ u_N\\ \end{pmatrix}. $$
All you now need is the vanishing of the determinant of the block matrix, all of whose blocks mutually commute, $$ \det \begin{pmatrix} -\lambda \mathbb{I} && \mathbb{I}\\ D && -\lambda \mathbb{I} \\ \end{pmatrix} = \det (\lambda ^2 \mathbb{I}- D)=0. $$ so the 2N eigenvalues of your (non-symmetric!) $2N\times 2N$ matrix are $\pm\sqrt{\lambda_1}$, ..., $\pm\sqrt{\lambda_N}$.
Leaving the eigenvectors to you... Check everything for N=2. Recall, for N =1, the eigenvectors $(1,\pm \sqrt \lambda_1 )^T$.