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.
1. 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.$
2. 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.$1 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. 1. 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.$2. 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\;\;$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.\$