This is only a partial answer, but it seems to be an interesting observation, and it's too long for the comments, so I write it here. The point of this answer is to show that the real difficulty lies in ensuring that the required condition holds for all $x$, because it is extremely easy to get a construction if the condition is only to be satisfied on a bounded interval, even if you require smoothness everywhere on $\mathbb R$ (and analyticity almost-everywhere). In fact, there is so little difficulty in making such a construction that any values of the $f(n)$ that you pick by hand will work (modulo some issues about convergence, so e.g. you can't pick values that diverge too quickly).
Claim. For any sufficiently nice sequence $(a_n)_{n\geq 1}$ (e.g. boundedness suffices), it is possible to find a smooth function
$f\colon\mathbb{R}\to\mathbb{R}$ satisfying $f(n)=a_n$ for positive
integers $n$, as well as $\displaystyle f(x)=\sum f(n)x^n$ in a
neighbourhood of zero.
To construct such a function we will exploit the fact that smooth functions are easy to bend around to our will, for example using bump functions. Let
$$ \sigma_{a,b}(x):=\begin{cases}1,&|x|\leq a\\\dfrac{\exp(-\frac{b-a}{x-a})}{\exp(-\frac{b-a}{x-a})+\exp(-\frac{b-a}{b-x})},&a<|x|<b,\\0,&|x|\geq b\end{cases}$$
be a smooth transition function from $1$ on the interval $(-a,a)$ to $0$ for $|x|\geq b$ for any real numbers $0<a<b$.
Proof. Let $f_r\colon (1-\epsilon,\infty)\to\mathbb R$ be any smooth function that interpolates the $a_n$, i.e. $f_r(n)=a_n$ for all $n\geq1$, for some fixed $\epsilon>0$. I claim that we can construct a function $f\colon\mathbb R\to\mathbb R$ which restricts to $f_r$ on $(1-\epsilon,\infty)$ with the desired properties. To do so consider the function on a neighbourhood of $0$ defined by $g\colon (-r,r)\to\mathbb R, g(x)=\sum a_nx^n$ (where $r$ is the radius of convergence of this series). Let $\epsilon'>0$ be any real number smaller than $\min\{r,1-\epsilon\}$. Then, the function $$f(x)=g(x)\sigma_{\epsilon',(1-\epsilon)}(x)+f_r(x)$$ has the desired properties with the identity $f(x)=\sum f(n)x^n$ holding in the neighbourhood $(-\epsilon',\epsilon')$. (Of course, the functions $g(x)$ and $f_r(x)$ are understood to mean $0$ outside their domains of definition.) $\square$
More generally, it is easy following this method to "twist the function to do whatever we want" in the range $(k,k+1)$ for any integer $k$. The difficulty is to "get past the integer values without ruining the desired property". Unfortunately, I have no idea how to do that.