Suppose on a manifold $M$ we have a differential operator $D$ of order 1 from the smooth sections $C^\infty(M, E)$ to smooth sections $C^\infty(M, F)$, where $E$ and $F$ are vector bundles of rank $n$ and $m$ respectively. If $u \in C^\infty(M, E)$ and $f \in C^\infty(M)$ is a smooth function, 
how would $D(fu)$ look like? I am trying to understand what the product/Leibniz rule could be in this case.



Generalized Leibniz rule
The "Generalized Leibniz rule" can be written in terms of commutators: see <http://en.wikipedia.org/wiki/Differential_calculus_over_commutative_algebras> and references therein.


In this perspective, the linear function in Stefan Waldmann's comment can be rewritten "without mentioning the $u$" as follows:
$$
\delta\_f(D):=[f,D]\in\textrm{Hom}\_{C^\infty(M)}(C^\infty(M,E),C^\infty(M,F))
$$
So, it is true that at this point you cannot say anything else.


Nevertheless, by "commuting again" you get zero:
$$
\delta\_g(\delta\_f(D))=0
$$
as "linearity" means "commuting with the multiplication". Now you can bring in the $u$ and get something resembling a Leibniz rule:
$$
D(fgu)= gD(fu)+fD(gu)-gfD(u).
$$
When $E$ is the trivial bundle, i.e., $C^\infty(M,E)=C^\infty(M)$, you can set $u=1$: then, $D$ is a derivation if and only if $D(1)=0$, and above formula become the usual Leibniz rule.


