A question regarding a cluster set

I've been having some difficulty trying to prove this exercise in $\textit{Complex Analysis}$ by Gamelin:

Let $\{z_k\}$ be a sequence of distinct points in a domain $D$ that accumulates at the boundary $\partial D$ and let $E$ be a closed subset of the extended complex plane $\mathbb{C} \cup \{\infty\}$. Show that there is an analytic function $f(z)$ on $D$ such that $E$ is the set of cluster values of $f$ along the sequence $\{z_k\}$.

Now there is a theorem that states that there is an analytic function $f(z)$ on the open unit disk $\mathbb{D}$ whose radial cluster set at any $\zeta \in \partial \mathbb{D}$ coincides with the entire complex plane. In the exercise couldn't I just let $D= \mathbb{D}$, and somehow work from there using this theorem? I'd appreciate some input here.

We can show rather more than what you ask for.

Let $z_k$ be a sequence of distinct points in an open subset $D$ of the extended complex plane $\mathbb{C}\cup\{\infty\}$, accumulating at the boundary $\partial D$. Let $a_k$ be a sequence in $\mathbb{C}$. Then, there is an analytic function $f$ on $D$ such that $f(z_k)=a_k$ for all $k$.

As any (nonempty!) closed subset of the extended complex plane is the set of cluster points of some sequence $a_k\in\mathbb{C}$, this gives what you ask.

To prove my claim, we can start by applying a Möbius transformation, if necessary, to transform the domain so that $\infty\in D$. Then, $\partial D$ is a compact subset of $\mathbb{C}$. Choose a sequence $y_k\in\partial D$ so that $\Vert z_k-y_k\Vert\to0$. The idea is to construct a function of the form $$ f(z)=\sum_{k=1}^\infty c_k(z-z_1)(z-z_2)\cdots(z-z_{k-1})\frac{(z_k-y_k)^{n_k}}{(z-y_k)^{n_k}} $$ where $c_k\in\mathbb{C}$ and $n_k\in\mathbb{N}$ are sequences which I will describe now. (Note, we will have $n_k\ge k-1$ whenever $c_k\not=0$, so that each of the terms in the sum is defined at $\infty$.)

Consider the partial sums, $$ f_m(z)=\sum_{k=1}^m c_k(z-z_1)(z-z_2)\cdots(z-z_{k-1})\frac{(z_k-y_k)^{n_k}}{(z-y_k)^{n_k}}. $$ If we have already chosen $c_k,n_k$ for $k < m$ such that $f_{m-1}(z_k)=a_k$ (all $k < m$), then it is automatically true that $f_m(z_k)=a_k$ (all $k < m$). There is also a unique choice for $c_m$ such that $f_m(z_m)=a_m$. Furthermore, the value of $c_m$ does not depend on how we choose $n_m$. Next, the term $$ c_m(z-z_1)(z-z_2)\cdots(z-z_{m-1})\frac{(z_m-y_m)^{n_m}}{(z-y_m)^{n_m}} $$ converges uniformly to zero on the set $S_m=\{z\in\mathbb{C}\colon\Vert z-y_m\Vert\ge2\Vert z_m-y_m\Vert\}$ in the limit $n_m\to\infty$. So, choosing $n_m$ large enough, we can suppose that this term is bounded by $2^{-m}$ on $S_m$.

We are now done. Every closed subset of $\mathbb{C}$ disjoint from $\partial D$ will be in $S_m$ for large $m$, so the sum defining $f$ converges uniformly on compacts in $D$ and $f$ is a well-defined analytic function. By construction, $f_m(z_k)=a_k$ for all $m\ge k$, so $f(z_k)=a_k$.