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}\prod_{j\neq i}((x-x_i)^n-(x_j-x_i)^n)^n\;\;$ $ 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.$
