We prove Runge type approximation results for linear partial differential operators with constant coefficients on spaces of smooth Whitney jets. Among others, we characterize when for a constant coefficient linear partial differential operator P(D) and for closed subsets \(F_1\subset F_2\) of \(\mathbb {R}^d\) the restrictions to \(F_1\) of smooth Whitney jets f on \(F_2\) satisfying \(P(D)f=0\) on \(F_2\) are dense in the space of smooth Whitney jets on \(F_1\) satisfying the same partial differential equation on \(F_1\) . For elliptic operators we give a geometric evaluation of this characterization. Additionally, for differential operators with a single characteristic direction, like parabolic operators, we give a sufficient geometric condition for the above density to hold. Under mild additional assumptions on \(\partial F_1\) and for \(F_2=\mathbb {R}^d\) this sufficient conditions is also necessary. As an application of our work, we characterize those open subsets \(\Omega\) of the complex plane satisfying \(\Omega =\operatorname {int}\overline{\Omega }\) for which the set of holomorphic polynomials are dense in \(A^\infty (\Omega )\) , under the additional hypothesis that \(\overline{\Omega }\) satisfies the strong regularity condition. Furthermore, for the wave operator in one spatial variable, a simple sufficient geometric condition on \(F_1, F_2\subset \mathbb {R}^2\) is given for the above density to hold. For the special case of \(F_2=\mathbb {R}^2\) this sufficient condition is also necessary under mild additional hypotheses on \(F_1\) .