I am looking for a reference for the following statement:

Assume that $P_1, \dots, P_k \in \mathbb R[x_1, \dots, x_m]$ and consider a system of PDEs \begin{align} P_i(\partial / \partial x_1, \dots, \partial / \partial x_m)u = f_i, \quad i = 1, \dots, k \end{align} where $u(x_1, \dots, x_m)$ is a scalar function to be found.

Then, if this system has a solution, we should obviously have $$ \sum_{i=1}^k Q_i(\partial / \partial x_1, \dots, \partial / \partial x_m)f_i = 0 $$ for each $k$-tuple $(Q_1, \dots, Q_k)$ such that $$ \sum_{i=1}^k Q_iP_i = 0. $$ The statement is that this condition is also sufficient for local solvability of our system of PDEs. More precisely,

1) if $f_1, \dots, f_k$ are $\mathrm C^\infty$ functions in a neighborhood of $x \in \mathbb R^m$, then there exists a $\mathrm C^\infty$ solution in a (possibly smaller) neighborhood of $x$;

2) if $f_1, \dots, f_k$ are analytic, then the solution is also analytic.

I suspect that this should be Ehrenpreis or Malgrange, but I was not able to find the precise statement.

Another question is whether there is a geometric way to understand this result. For instance, if $P_1, \dots, P_k$ are of degree $1$, then the statement can be deduced from the Frobenius theorem.

there existsa solution $u$ that is analytic'. (Obvious examples show that notallof the solutions need be analytic, even when the $f_i$ are all zero.) Such a version of the statement (at least for local solutions) follows from the Cartan-Kähler Theorem. $\endgroup$