The problem may be split into two independent and classical ones: the Hermite interpolation, and the Weierstrass approximation.

First, we want a polynomial $p\in \mathbb{R}[x]$ with given derivatives at some given nodes $x _ 0,\dots, x _ m $. This is an instance of the Hermite interpolation problem; yours has exactly one solution $p$ with $ \operatorname{deg}(p) < (m+1)(n+1) $ (the degree one would expect in terms of number of linear conditions).
So, given your $f\in C^n$, you can find a polynomial with $p^{(j)}(x _ i)=f^{(j)}(x _ i)$ for all $0 \le i \le m$ and $0 \le j \le n$.

Second, as a consequence,

$$\frac{f(x)-p(x)}{\prod _ {i=0}^ m(x- x _ i)^n}$$
is (extends to) a continuous function on $[a,b]$ that vanishes on the points $x _ i$. By the Stone-Weierstrass approximation theorem there is a polynomial vanishing on the points $x _ i$ as well, whose uniform distance from that function on the interval $[a,b]$ is less than, say, $\epsilon (b-a)^{-n(m+1)}$. In other words, there is a polynomial $q\in \mathbb{R}[x]$ such that

$$\bigg\| \frac{f(x)-p(x)}{\prod _ {i=0}^ m(x- x _ i)^n} - q(x) \prod _ {i=0} ^ m(x- x _ i) \bigg\|_{\infty, [a,b]} < \epsilon (b-a)^{-n(m+1)}\ , $$

therefore the polynomial $P(x):= p(x)+ q(x) \prod _ {i=0} ^ m(x- x _ i) ^{n+1}$ fullfills the requirements, for $P^{(j)}(x _ i)=p^{(j)}(x _ i)=f^{(j)}(x _ i)$ for all $0 \le i \le m$ and $0 \le j \le n$, and
$$\|f- P\|_{\infty,[a,b]} < \epsilon\ .$$

**btw.** Incidentally, some time ago I happen to notice that one can find the solution of the Hermite interpolation problem as an application of the Chinese Remainder Theorem in the ring of polynomials, and wrote here the details.

**edit.** As to why The set $A$ of all polynomial functions on $[a,b]$ that vanish at given points $x_0,\dots, x_m$ is dense in all continuous functions on $[a,b]$ that vanish in $x_0,\dots, x_r$. One way, a bit abstract but quite immediate is, to see it as a corollary of the Stone-Weierstrass theorem (A separating closed algebra of real valued functions on a compact space $X$ is either $C(X)$ or a maximal ideal $M_x\subset C(X)$, the set of all functions vanishing at $x$). Consider $X=$ the topological quotient of $[a,b]$ obtained identifying all points $x_i$ to a point $\xi$. All functions in $A$ factor through to the quotient map, and define a closed separating algebra of continuous functions on $X$ that vanish on the identified point $\xi$. Thus, this algebra contains all continuous functions on $X$ that vanish on $\xi$, which is the thesis read on the quotient.

Note that the same construction holds in general, and provides a characterization of all closed algebras $A$ of continuous functions on a compact space $X$: identifying all points that are not distinguished by the functions of $A$ (that is, under the equivalence relation $x R_A y$ iff $f(x)=f(y)$ for all $f\in A$) one gets a Hausdorff compact quotient space (whether or not $X$ is Hausdorff), and the quotient map $\pi: X\to X/{R _ A}$ induces an isometric isomorphism of algebras $f\mapsto f\circ \pi$ of either $C(X/{R_A})$ or a maximal ideal of it onto $A$; conversely, any Hausdorff quotient of $X$ produces a closed sub-algebra of $C(X)$ this way.

Another way to see it is as a corollary of the classic Weierstrass theorem: Consider $P$ as in your comment below, then add a perturbation $L$ that makes $P+L$ vanish on the points $\{x_i\}$; this has been clearly explained in Ilya Bogdanov's answer. Here you don't have derivatives and $L$ is just a Lagrange interpolation polynomial, which is small in the uniform norm because it is small on the points $\{x_i\}$.