Since Giuseppe has already given an answer, I will give the details of the extension.
Let $A\subset \mathbb R$ a set which consist of isolated points and for each $a\in A$, $\{x(a)_n\}$ a sequence of real numbers. We can find a smooth function $f$ such that for all $a\in A$ and $n\in\mathbb N$, we have $f^{(n)}(a)=x(a)_n$.
Since $\mathbb R$ is separable, this set has to be countable, hence we can write $A=\{a_n,n\in\mathbb N\}$. We can write $b_{n,k}=x(a_n)_k$ for each $n$ and $k$. For each $n$, we consider $r_n $such that $\left]a_n-3r_n,a_n+3r_n\right[\cap A=\{r_n\}$. Let $g_n$ be a bump function, $g_n=1$ on $\left]-r_n,r_n\right[$ and $\mathrm{supp}g_n\subset \left[-2r_n,2r_n\right]$. Put $f_{n,j}(x)= \frac{b_{n,j}}{j!}(x-a_n)^jg_n\left(\frac{x-a_n}{\varepsilon_{n,k}}\right)$ where $\varepsilon_{n,k}=\frac 1{1+|b_{n,k}|}\frac 1{4^kM_{n,k}}$, with $\displaystyle M_{n,k}:=\max_{0\leq j\leq k}\sup_{x\in\mathbb R}|g_n^{(j)}(x)|$. Note that $f_{n,k}$ has a compact support, and putting for all $n$: $\displaystyle h_n(x):=\sum_{k=0}^{+\infty}f_{n,k}(x)$, $h_n$ is a smooth function. Indeed, for $n$ fixed and $k\geq d+1$ we have $\sup_{x\in\mathbb R}|f_{n,k}^{(d)}|\leq \frac 1{2^k}$, hence the normal convergence of the series on the whole real line for all the derivatives. A (boring but not hard) computation which uses Leibniz rule gives that we have $h_n^{(k)}(a_n)=b_{n,k}$.
Now, put $f(x):=\sum_{n=0}^{+\infty}\sum_{j=0}^{+\infty}f_{n,j}(x).$ Since the supports of $h_n$ are pairwise disjoint, the series in $n$ is convergent to a smooth function.
We will note that we cannot extend the result to sets $S$ which contain a point which is not isolated, because the continuity of the derivative gives a restriction on the sequence $\{x(s)_n\}$ for the $s$ in a neighborhood of the non-isolated point. Namely, if $s_0$ is a non-isolated point and $\{s_n\}\subset S$ a sequence which converges to $s_0$ then the sequence $\{x(s_n)_0\}$ has to converge to $x(s_0)_0$ (hence we can't get all the sequences).