Crank-Nicolson for coupled PDE's $\newcommand{\T}{T}$
$\newcommand{\partiald}[2]{\frac{\partial #1}{\partial #2}}$
$\newcommand{\partialdd}[2]{\frac{\partial^2 #1}{\partial #2^2}}$
I am trying to solve a set of coupled PDE's with the Crank-Nicolson method. So far I have used it to solve a single PDE, the 1D diffusion problem in the Wikipedia article I have linked. This involved turning the set of equations into a matrix equation
$$
\begin{bmatrix}
A & 1 & 0 & \cdots &     & 0      \\
1 & A & 1 & \cdots &     &        \\
0 & 1 & A & \ddots &     &        \\
\vdots & \vdots & \ddots & \ddots \\
  &   &   &        &  A  & 1      \\
0 &   &   &        &  1  & A
\end{bmatrix}
\begin{bmatrix}
\T^{n+1}_1\\
\T^{n+1}_2\\
\T^{n+1}_3\\
\vdots \\
\T^{n+1}_{N-2}\\
\T^{n+1}_{N-1}
\end{bmatrix}
=
\begin{bmatrix}
d_1^n - \T_0^{n+1}\\ 
d_2^n \\
d_3^n \\
\vdots \\ 
d_{N-2}^n\\
d_{N-1}^n - \T_{N}^{n+1}
\end{bmatrix} 
$$
where $A = -\frac{1+2r}{r}$ and 
$$d_i^n = - \left[\T_{i+1}^{n} + \left(\frac{1-2r}{r}\right) \T_i^{n} + \T_{i-1}^{n}\right]$$
I specify $T_0^{n+1}$ and $T_N^{n+1}$ as boundary conditions, and provide an initial set of values for $T_i^0$. From there it is simply a matter of iteratively solving the set of linear equations for the $T_i^{n+1}$.
Now I want to apply this method to solve the set of coupled equations
\begin{align*}
 \partiald{v}{t} &= - \partiald{\rho}{x} \\
 \partiald{\rho}{t} &= - \partiald{v}{x} 
\end{align*}
Applying the Crank-Nicolson algorithm I find the following relation
\begin{align*}
  v_i^{n+1} + r\left(\rho_{i+1}^{n+1} - \rho_{i-1}^{n+1}\right) &= v_i^n - r\left(\rho_{i+1}^n - \rho_{i-1}^n\right) \\
  \rho_i^{n+1} + r\left(v_{i+1}^{n+1} - v_{i-1}^{n+1}\right) &= \rho_i^n - r\left(v_{i+1}^n - v_{i-1}^n\right)
\end{align*}
where $r \equiv \Delta t/ 4\Delta x$. But from here I am stuck, I don't know how to turn this into a matrix equation due to having two variables. I am assuming I want to turn this into a vector equation, like the way I did in this answer, but I'm not sure how in this case I would then convert it into a matrix equation which can be iteratively solved.
Edit:
I realize that I can convert these into two separte equations
\begin{align*}
\partialdd{\rho}{x} &= \partialdd{\rho}{t}\\
\partialdd{v}{x} &= \partialdd{v}{t}
\end{align*}
which I think I can solve using the previous method, though I think I will need to use an explicit algorithm for the first time step. But my question above still stands as to how to solve the set of coupled equations.
 A: $\renewcommand{\d}{\vec{d}}
\renewcommand{\S}{\vec{S}}
\renewcommand{\v}{v}
\renewcommand{\r}{\rho}$
I solved the problem by doing the following: starting from
\begin{align*}
  \v_i^{n+1} + r\left(\r_{i+1}^{n+1} - \r_{i-1}^{n+1}\right) &= \v_i^n - r\left(\r_{i+1}^n - \r_{i-1}^n\right) \\
  \r_i^{n+1} + r\left(\v_{i+1}^{n+1} - \v_{i-1}^{n+1}\right) &= \r_i^n - r\left(\v_{i+1}^n - \v_{i-1}^n\right)
\end{align*}
Define a vector $\vec{S}_i^n = (\tilde{v}_i^n , \r_i^n)^{\mathrm T}$. Then we can write these two equations as
$$
\S_i^{n+1}
+
\begin{bmatrix}
0 & r \\ r & 0
\end{bmatrix}
\left(\S_{i+1}^{n+1} - \S_{i-1}^{n+1}\right)
= 
\S_i^n 
-
\begin{bmatrix}
0 & r \\ r & 0
\end{bmatrix}
\left(\S_{i+1}^n - \S_{i-1}^n\right)
$$
call this matrix $A$. Then we write our set of equations as
$$
\begin{bmatrix}
1  & A      &        &       &  0  \\
-A & 1      & \ddots &       &     \\
   & \ddots & \ddots &       &     \\
   &        &        &     1 &  A  \\
0  &        &        &    -A &  1  
\end{bmatrix}
\begin{bmatrix}
\S_1^{n+1}     \\
\S_2^{n+1}     \\
\vdots         \\
\S_{N-2}^{n+1} \\
\S_{N-1}^{n+1} 
\end{bmatrix}
=
\begin{bmatrix}
1 & -A     &        & \!   &  0  \\
A & 1      & \ddots & \!   &     \\
  & \ddots & \ddots & \!   &     \\
  &        &        & \! 1 & -A  \\
0 &        &        & \! A &  1  \\
\end{bmatrix}
\begin{bmatrix}
\S_1^{n}     \\
\S_2^{n}     \\
\vdots       \\
\S_{N-2}^{n} \\
\S_{N-1}^{n} \\
\end{bmatrix}
+
\begin{bmatrix}
A\left(\S_0^{n} + \S_0^{n+1}\right)  \\
0            \\ 
\vdots       \\
0            \\
-A\left(\S_N^{n} + \S_N^{n+1}\right)  \\
\end{bmatrix}
$$
