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$.
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) such that $$ 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 $$ 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.