Maybe I am asking a triviality. If this 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 existence of a solution $f$? (we can shrink $U$ if needed).

I don't think I can apply the [Cauchy–Kovalevskaya theorem](https://en.wikipedia.org/wiki/Cauchy%E2%80%93Kowalevski_theorem) in this case, because we have several vector fields...

On the other hand, at the end of [Wikipedia page of Cauchy–Kovalevskaya theorem](https://en.wikipedia.org/wiki/Cauchy%E2%80%93Kowalevski_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.