System of linear first order PDE with constant coefficients recently in my researches I've come across the following operator
$$L\left(\begin{array}{c}
         a_1\\
         \vdots\\
         a_n
         \end{array}\right)=M_1\left(\begin{array}{c}
         \partial_{z_1}a_1\\
         \vdots\\
         \partial_{z_1}a_n
         \end{array}\right)+\dotso+M_n\left(\begin{array}{c}
         \partial_{z_n}a_1\\
         \vdots\\
         \partial_{z_n}a_n
         \end{array}\right),$$
where $a_1,\dotso,a_n:\mathbb R^n\to\mathbb R$ are smooth functions of the independent variables $z_1,\dotso,z_n$ and $M_1,\dotso M_n$ are real square $n\times n$ matrices with constant coefficients (though none of them is invertible).
Is there a general theory to treat this simple system of linear first order pdes? Is it possible to write explicitly a solution if we were to solve
$$L\left(\begin{array}{c}
         a_1\\
         \vdots\\
         a_n
         \end{array}\right)=\left(\begin{array}{c}
         b_1\\
         \vdots\\
         b_n
         \end{array}\right),$$
for some smooth functions $b_1,\dotso,b_n:\mathbb R^n\to\mathbb R$?
Thank you for the references and the patience.
Kind regards,
Guido
Edit
Following the suggestions, I report here the $3\times 3$ system which is the toy model I'm currently studying
$$\left(\begin{array}{ccc}
                 B_{11} & B_{12} & 0\\
                    -B_{11} & 0 & 0\\
                 0 & -B_{12} & 0
              \end{array}
           \right)\left(\begin{array}{c}
\partial_{z_1}a_1\\
\partial_{z_1}a_2\\
\partial_{z_1}a_3\\
                           \end{array}\right)+
\left(\begin{array}{ccc}
                 0 & 0 & -B_{23}\\
                    0 & B_{22} & B_{23}\\
                 0 & -B_{22} & 0
              \end{array}
           \right)\left(\begin{array}{c}
\partial_{z_2}a_1\\
\partial_{z_2}a_2\\
\partial_{z_2}a_3\\
                           \end{array}\right)+
\left(\begin{array}{ccc}
                 0 & 0 & -B_{33}\\
                    -B_{31} & 0 & 0\\
                 B_{31} & 0 & B_{33}
              \end{array}
           \right)\left(\begin{array}{c}
\partial_{z_3}a_1\\
\partial_{z_3}a_2\\
\partial_{z_3}a_3\\
                           \end{array}\right)=
\left(\begin{array}{c}
      b_1\\
      b_2\\
      -b_1-b_2
      \end{array}\right),
$$
$a_i$ and $b_i$ are smooth functions, the $b_i$'s are compactly supported. I'm interested in solving locally the above system. I think it has something to do with the usual theory for the hyperbolic systems of differential equations, however I'm to unexperienced in the field to be sure about the steps to follow. I understand the condition to be imposed on the functions $b_i$ as a necessary condition for the system to be possibly solvable, but it is not clear to me why there should exist a solution at all! This is the most I can say up to now. Thanks again for the patience and the help
 A: Since you wrote the system explicitly, here's an explicit solution using Fourier transforms. Let $a(z) = \int e^{ik\cdot z} a(k)\, d^3k $ and $b(z) = \int e^{ik\cdot z} b(k)\, d^3k$. The $a(k)$ and $b(k)$ vectors will satisfy the Fourier transformed equation $iM(k) a(k) = b(k)$, where $M(k) = (M_1 k_1 + M_2 k_2 + M_3 k_3)$. It can be checked explicitly that the following formula gives a complete solution for $a(k)$:
$$\begin{gathered}
\begin{pmatrix} a_1(k) \\ a_2(k) \\ a_3(k) \end{pmatrix} =
\frac{i}{\Delta(k)}
\begin{pmatrix}
  B_{23} k_2 & B_{23} k_2 + B_{33} k_3 \\
  0 & 0 \\
  B_{11} k_1 + B_{31} k_3 & B_{11} k_1
\end{pmatrix}
\begin{pmatrix} b_1(k) \\ b_2(k) \end{pmatrix}
+ \begin{pmatrix} \Delta_1(k) \\ \Delta_2(k) \\ \Delta_3(k) \end{pmatrix} c(k), \\
\Delta(k) = B_{11} B_{33} k_1 k_3 + B_{23} B_{31} k_2 k_3 + B_{31} B_{33} k_3^2 , \\
\Delta_1(k) = B_{12} B_{23} k_1 k_2 + B_{22} B_{23} k_2^2 + B_{22} B_{33} k_2 k_3 , \\
\Delta_2(k) = B_{11} B_{33} k_1 k_3 + B_{23} B_{31} k_2 k_3 + B_{31} B_{33} k_3^2 , \\
\Delta_3(k) = B_{11} B_{12} k_1^2 + B_{11} B_{22} k_1 k_2 + B_{12} B_{31} k_1 k_3 ,
\end{gathered}
$$
where $c(k)$ is an arbitrary scalar function.
Of course, because the denominator $\Delta(k)$ vanishes for some $k$, the Fourier coefficients $a(k)$ must be interpreted as distributions.
A: Perhaps this is too naive an answer, but since this problem has constant coefficients, it seems perfectly set up for separation of variables. Look for complex solutions of the form $e^{i\vec{\lambda}\cdot\vec{x}} \vec{A}(\vec{\lambda})$. The first equation then takes the form $0=(-L+\sum M_i \lambda_i)\vec{A}$. This has a nonzero solution iff $det(-L+\sum M_i\lambda_i)=0$. The determinant is a polynomial in $\vec{\lambda}$, so there is a well established theory for finding the solution surfaces. (Of course, in high dimensions, well established does not necessarily mean easy.) Once the frequencies $\vec{\lambda}$ have been found, it is an eigenvalue problem to find the associated solutions $\vec{A}$. Since the original matrices were real valued, the frequencies $\vec{\lambda}$ will either be real or come in complex pairs, allowing for the construction of real solutions. 
Be aware that series solutions of this type conceal a lot of information.  By appropriate choices of $M_i$, each component of $\vec{a}$ can be made to satisfy the Laplace equation or to satisfy the wave equation. These have radically different behaviours. The Laplace equation has no nonzero solutions that vanish at infinity, whereas the wave equation has an infinite dimensional space of such solutions (for $n>1+1$). 
As the comment with curl above shows, introducing $\vec{b}=L\vec{a}$ does not seem to help the analysis. 
