I have a linear, first-order homogeneous PDE system with polynomial coefficients $$L_j\, f =0,\text{ for } j=1,..,J\quad \text{ where } L_j \text{ is a first order, diff. operator with polynomial coeff.}$$ so any constant function always solves it. I want to find whether in a neighborhood of a point p, all solutions are constant or not.
We compute the Lie-brackets, $L_{ij}=[L_i,L_j]$, and since $L_{ij}\,f =0$ provided $f$ is a solution of the original system, I can think of $L_{ij}$ as a vector field in the tangent bundle of $f$. If we get the minimal module of linear operators that contains the $L_i$ and is "closed" under Lie-bracket, evaluate them at p, and the compute the dimension of the corresponding vector space (which is a subspace of the tangent space of a solution at p) we have that:
All solutions are constant at p if and only if the space described above has maximal dimension. Is this true? I think this statement is correct, if not please let me know.
I think the tangent bundle you should consider is spanned by the Lie algebra generated by $L_i$, not just $L_{ij}$. In particular, if the Lie algebra is commutative then $L_{ij}=0$, which is not what you are looking for.
My favourite reference on these things is "Mathematical Methods of Classical Mechanics" by V.I. Arnold. In the chapter on Hamiltonian Mechanics he discusses the differential operators that preserve the flow and their Poisson bracket (which makes them into Lie algebra, as you described). Each $L_i$ determines a "first integral" of the flow, and if you can find $k$ of them in involution (that is, $L_{ij}=0$) then you have $k$ invariant quantities. You claim is that if in $n$-dimensional phase space you have $k=n$ then the solutions are constant, and this is true. $k=n$ is very rare though; it's more common to have $k<n/2$. However, even for $k=n/2$ there is a very good result: with $k=n/2$ first integrals in involution the system is completely solvable.