Here is a somewhat brute-force, but direct construction. We must require that the $a_n$s are all different, of course. Without loss of generality we can also assume that $a_n\ne 0$ for $n>0$.
We will construct a sequence of entire functions that have the right value for increasing prefixes of the $a_n$s, and then show that the functions converge to an entire function.
In step $0$, set $f_0(z)=b_0$.
In step $n$ for $n>1$, we construct $f_n$ as $f_{n-1}+g_n$, where $g_n$ is a function that is $0$ at $a_0,\ldots,a_{n-1}$ and has the right value at $a_n$ to make $f_n(a_n)=b_n$. In particular we use
$$ g_n(z)=cz^k\prod_{i=0}^{n-1}(z-a_n) $$
for an appropriate constant $c$ and a $k$ chosen large enough that $\lvert g(z)\rvert \le 2^{-n} $ whenever $\lvert z\rvert \le \lvert a_n\rvert-1$.
Now, for any $D>0$, $\lvert a_n\rvert-1$ will eventually be greater in $D$, and therefore $ \lim_n f_n = \sum_{n=0}^\infty g_n $ will converge uniformly for $\lvert z\rvert < D$. But all of the $f_n$ are entire (in fact polynomials), and uniform convergence on an open subset of $\mathbb C$ preserves analyticity, so the limit is analytic on $\lvert z\rvert < D$. But $D$ was arbitrary, so $\lim_n f_n$ is analytic everywhere.
