Existence of solution to linear inhomogeneous first order PDEs systems Maybe I am asking a triviality. If that is the case, please let me know and I will close the question. I have searched a lot but I didn't find an adequate and forceful response.
For $i=1,\ldots, r$, let $Z_i$ be $r$ linearly independent vector fields defined on an open subset $U$ of $\mathbb{R}^n$, $r<n$. And let $\lambda_i$ be $r$ smooth functions defined on $U$. We can consider the inhomogeneous system of first-order linear PDEs
$$
Z_i(f)=\lambda_i, \quad i=1,\ldots,r.
$$
Question: What hypothesis do we need to assure the local existence of a solution $f$?
I don't think I can apply the Cauchy–Kovalevskaya theorem in this case, because we have here several vector fields...
On the other hand, at the end of Wikipedia page of Cauchy–Kovalevskaya theorem it is named the Cauchy–Kovalevskaya–Kashiwara theorem, which seems very technical to me, but I feel it has something to do...
Finally, I think that involutivity of the distribution $\{Z_1,\ldots,Z_r\}$ plays an important role here. For example, in case $\lambda_i=0$ Frobenius theorem would tell us that there exist $n-r$ functionally independent solutions.
 A: You are correct that Cauchy-Kovalevskaya does not apply directly to this problem, but there are other theorems that give sufficient conditions, provided that you make certain basic regularity assumptions.
For example, in the smooth involutive case, i.e., when the $Z_i$ (as well as the $\lambda_i$) are also sufficiently differentiable and satisfy the involutivity condition that $[Z_i,Z_j] = c^k_{ij}\,Z_k$ (summation convention assumed) for some smooth functions $c^k_{ij}=-c^k_{ji}$, then an obvious necessary condition on the $\lambda_i$ for local solvability of the system $Z_if = \lambda_i$ is
$$
Z_i(\lambda_j)-Z_j(\lambda_i) = c^k_{ij}\lambda_k\,\qquad \forall i,j
\quad\,1\le i,j\le r,\tag1
$$
since both sides would equal $[Z_i,Z_j]f$ for any (local) solution $f$.  It's a consequence of the Pfaff-Darboux Theorem (appropriately interpreted) that, if condition (1) is satisfied, then, every $p\in U$ has an open neighborhood $V\subset U$ on which there exists an $f$ satisfying the given system $Z_if = \lambda_i$.  However, because the system is linear inhomogeneous, there is a simpler proof of this result that only uses the Frobenius Theorem.  I give that proof in a remark at the end of this answer.
This can be extended to more general situations:  Suppose that the distribution $D$ spanned by the $Z_i$ is not involutive, but that we can choose new $Z_a$ for $r<a\le r'$ such that $Z_1,\ldots, Z_{r'}$ are linearly independent and give a basis for the sections of the distribution $D'$ of constant rank $r'$ spanned by the $Z_i$ and $[Z_i,Z_j]$ for $1\le i,j\le r$.  Then writing $[Z_i,Z_j] = c^a_{ij} Z_a$ (where now, the index $a$ runs from $1$ to $r'$), we see that a necessary condition for solvability is that there exist functions $\lambda_a$ for $r<a\le r'$ (necessarily unique) that satisfy
$$
Z_i(\lambda_j)-Z_j(\lambda_i) = c^a_{ij}\lambda_a\,\qquad \forall i,j\qquad\,1\le i,j\le r,\tag2
$$
(where the implied summation on $a$ runs from $a=1$ to $a=r'$)
and that we would have to have $Z_af = \lambda_a$ for $1\le a\le r'$.
This is now a bigger system (the 'prolongation' of the original system), and we can check whether $D'$ is involutive.  If not, we expand (aka,
'prolong') the system again.  If we can repeat this prolongation process until we reach a system $D''$ for some $r''>r$ that is involutive, then we can apply the first criterion to determine local solvability of the original system.
Along the way, we might run into situations where the new $\lambda_a$ one needs to find don't exist on some open set, in which case, there won't be a solution to the original system.  More complicated things can happen if the 'prolonged' distributions don't have constant rank, but, for most practical purposes, the above process will give an effective test as to whether and where in $U$ (local) solutions $f$ exist.
Remark:  I can sketch the proof of the basic existence theorem mentioned above, based on the Frobenius theorem.
First, suppose that $[Z_i,Z_j] = 0$ for $1\le i,j,\le r$.  Then, by the Simultaneous Flow Box Theorem, for each $p\in U$ there exists an open $p$-neighborhood $V\subset U$ on which there exists a $p$-centered coordinate chart $z = (z^1,\ldots,z^n)$ such that $Z_i = \partial/\partial z^i$ on $V$ for $1\le i\le r$. On $V$, the given system becomes $\partial f/\partial z^i = \lambda_i$, and the condition (1) becomes $\partial\lambda_i/\partial z^j = \partial\lambda_j/\partial z^i$, and, by the usual integration formula, it is clear that, when this holds, there will be a smaller $p$-neighborhood $V'\subset V$, on which there exists a function $f$ satisfying $\partial f/\partial z^i = \lambda_i$ for $1\le i\le r$.
Now consider the more general case in which $[Z_i,Z_j] = c^k_{ij}\,Z_k$.  By the Frobenius Theorem, every point $p\in U$ will have a $p$-neighborhood $V$ with a $p$-centered coordinate system $z = (z^1,\ldots,z^n)$ such that $Z_i = a_i^k\,\partial/\partial z^k$ on $V$ for some smooth functions $a_i^k$ on $V$. Moreover, the $r$-by-$r$ matrix $a = (a_i^k)$ is invertible because the $Z_i$ for $1\le i\le r$ are linearly indepdendent on $U$.  Consequently, we can write $\lambda_i = a_i^k\mu_k$ for some functions $\mu_k$ ($1\le k\le r$) on $V$.
Now, we can expand $[Z_i,Z_j] = c^k_{ij}\,Z_k$ to get a formula for $c^k_{ij}$ in terms of the $a_i^k$ and their partials with respect to the $\partial/\partial z^j$.  Using this formula and some index juggling, we see that
$$
Z_i(\lambda_j)-Z_j(\lambda_i) - c^k_{ij}\lambda_k
= a_i^ka_j^l\left(\frac{\partial \mu_k }{\partial z^l}-\frac{\partial \mu_l}{\partial z^k}\right).
$$
Since the given system $Z_if=\lambda_i$ is equivalent to $\partial f/\partial z^k = \mu_k$, we see that we are reduced to the case $c^k_{ij} = 0$, which has already been treated.
