Approximation by polynomials Let $f:[a,b] \rightarrow \mathbb{R}$ be of class $C^n$. Let $ x_0, ..., x_m$ be different numbers from $[a,b]$. 
Does for each $\varepsilon >0$ there exist a polynom $P$ such that $P^{(k)}(x_i)=f^{(k)}(x_i)$ for $i=0,...,m$, $k=0,...,n$ 
and $sup_{x \in [a,b]} |f(x)-P(x)|< \varepsilon$?
 A: For the sake of simplicity, let us assume that $[a,b]=[0,1]$. $C^k$ always means $C^k[0,1]$. We will even approximate $f$ in $C^n$-norm satisfying your additional condition.


*

*As was mentioned in the comments, you can easily approximate $f$ together with all its derivatives up to $n$th uniformly by a polynomial. In fact, it is enough to approximate $f^{(n)}$ with an adequate accuracy: if $||f'-P'||_C<\varepsilon$ and $f(0)=P(0),$ then $||f-P||_C<\varepsilon.$

*Now take the polynomials $Q_{ik}(x)$ such that $Q_{ik}^{(d)}(x_j)=0$ for all $d=0,\dots,n$ 
and $j=0,\dots,m$ except that $Q_{ik}^{(k)}(x_i)=1$. Such polynomials are easy to construct: for instance, one may take 
$$
  Q_{ik}(x)=c_{ik}(x-x_i)^k\prod_{j\neq i}\left((x-x_i)^{n+1}-(x_j-x_i)^{n+1}\right)^{n+1}\;\;
$$ 
for a suitable constant $c_{ik}.$ Let $M=\max_{i,k}||Q_{ik}||_{C^n}.$ Then, let the approximation in the previous paragraph be $\delta$-accurate with $\delta=\varepsilon/(2M(m+1)(n+1)).$ To correct the values of the polynomial and its derivatives at $x_i,$ it is enough to add the polynomials $Q_{ik}$ multiplied by the coefficients with absolute values $\leq\delta,$ hence the total error will be not more that $\delta+(m+1)(n+1)M\delta<\varepsilon.$
A: 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\}$. 
