On finding polynomials that approximate a function and its derivative (extensions of Stone-Weierstrass?) The Stone-Weierstrass Theorem tells us that we can approximate any continuous $f:\mathbb{R}^n\to\mathbb{R}$ arbitrary well on a compact subset of $\mathbb{R}^n$ by some polynomial. Suppose that $f$ is continuously differentiable, under what conditions can we guarantee we can find a polynomial $p$ such that $p$ approximates arbitrarily well $f$ and the partial derivatives approximate of $p$ arbitrarily well those of $f$ (on a compact subset of $\mathbb{R}^n$)? References are welcome. See below for my (probably naive) approach to the question.

It is straightforward to find conditions if $n=1$: Let $K:=[a,b]\subseteq\mathbb{R}$ and let $K_+$ denote some open cover of $K$. Suppose that $f:K_+\to\mathbb{R}$ is continuously differentiable and let $f'$ denote its derivative. Then, by Weierstrass's approximation Theorem we have that there exists a polynomial $q$ such that 
$$||f'(x)-q(x)||\leq \varepsilon$$
for all $x\in K$. Choose some $y\in K$. For any $x\in K$ define
$$p(x):=f(y)+\int_y^xq(t)dt.\quad\quad(*)$$
By the Fundamental Theorem of calculus we have that $p'=q$. So for any $x\in K$
$$||f(x)-p(x)||=\left|\left|\int_y^x (f'(t)-q(t))dt\right|\right|\leq \int_y^x ||f'(t)-q(t)||dt\leq \varepsilon||x-y||\leq \varepsilon (a-b).$$
However, the above argument doesn't work in higher dimensions because the analogue of $(*)$ does not imply that $\nabla p = q$ (since not every vector of polynomials spans a gradient field). Can anyone suggest a way around this? That is, a way to to show that for continuously differentiable $f:\mathbb{R}^n\to\mathbb{R}$ and for any $\varepsilon>0$, we can always find a polynomial $p$ such that 
$$||f(x)-p(x)||+||\nabla f(x)-\nabla p(x)||\leq \varepsilon$$
for all $x$ in some compact set $K$? Feel free to add extra conditions on $K$, if it makes life easier.
 A: By using a partition of unity, we may assume that $f$ has compact support. Let $\phi_m$ be the density of the normal vector with mean $\bf 0$ and covariance matrix $I/m$. Let $q_m$ be a Maclaurin polynomial approximating $\phi_m$ of sufficiently high degree that $\lim_{m \to \infty} \|\phi_m - q_m\|_{C(K)} = 0$, where $K$ is the Minkowski difference of where you want the approximation to hold and the support of $f$. Now let $p_m := q_m * f$. These $p_m$ do the job.
Remark: Nachbin has a generalization of the Stone-Weierstrass theorem to include approximating derivatives, even on manifolds. See MR0030590: 
Nachbin, Leopoldo, Sur les algèbres denses de fonctions différentiables sur une variété. C. R. Acad. Sci. Paris 228, (1949), 1549–1551. 
A: If we assume a little more smoothness, we can use a trick similar to the one in this question.
Let me take $n=2$ for concreteness.  Suppose $f$ is $C^2$ on a neighborhood $U$ of $K$.  With an appropriate cutoff function we can find a $C^2$, compactly supported $g$ which agrees with $f$ on a (possibly smaller) neighborhood of $K$.  We can also find a large enough square $Q = (-M,M)$ which contains the support of $g$, so that $g$ and all its derivatives vanish on the boundary $\partial Q$.
Fix $\epsilon > 0$.  The second partial $\partial_x \partial_y g$ is continuous, so we can find a polynomial $p_0$ with $|p_0 - \partial_x \partial_y g| < \epsilon$ on $\bar{Q}$.  Set
$$p(x,y) = \int_{-M}^x \int_{-M}^y p_0(s,t)\,dt\,ds.$$
Then by the fundamental theorem of calculus and Fubini's theorem we have, for any $(x,y) \in Q$,
$$\begin{align*}|p(x,y) - g(x,y)| &= \left|\int_{-M}^x \int_{-M}^y p_0(s,t) - \partial_x \partial_y g(s,t)\,ds\,dt\right|\\ & \le  \int_{-M}^x \int_{-M}^y |p_0(s,t) - \partial_x \partial_y g(s,t)|\,ds\,dt \\&\le (2M)^2 \epsilon \end{align*}$$
and 
$$\begin{align*}|\partial_x p(x,y) - \partial_x g(x,y)| &= \left| \int_{-M}^y p_0(s,t) - \partial_x \partial_y g(s,t)\,dt\right|\\ & \le \int_{-M}^y |p_0(s,t) - \partial_x \partial_y g(s,t)|\,dt \\&\le 2M \epsilon. \end{align*}$$
A similar argument works for $\partial_y$.  Since $f=g$ on a neighborhood of $K$, we have that $p$ and $f$, as well as their derivatives, are uniformly close on a neighborhood of $K$.
In $\mathbb{R}^n$, this approach requires us to assume that $f$ is $C^n$.
