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 in this case, because we have 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.