I will assume we can write
$$f(z) = \frac{c}{z_0-z} + \sum_{n=0}^{\infty} b_n z^n$$
for some value of $c \ne 0$, and $\lim_{n \to \infty} b_n = 0$. Then
$$f(z) = \sum_{n=0}^{\infty} a_n z^n$$
where
$$a_n = b_n + \frac{c}{z_0^{n+1}}$$
Then
$$\begin{align}\lim_{n \to \infty} \frac{a_n}{a_{n+1}} &= \lim_{n \to \infty} \frac{\displaystyle b_n + \frac{c}{z_0^{n+1}}}{\displaystyle b_{n+1}+ \frac{c}{z_0^{n+2}}}\\ &= \lim_{n \to \infty} \frac{\displaystyle \frac{c}{z_0^{n+1}}}{\displaystyle \frac{c}{z_0^{n+2}}}\\ &= z_0\end{align}$$
as was to be shown. Note that the second step above is valid because $z_0$ is on the unit circle.
For a nonsimple pole, we may write
$$f(z) = \frac{c}{(z_0-z)^m} + \sum_{n=0}^{\infty} b_n z^n$$
for $m \in \mathbb{N}$. It might be known that
$$(1-w)^{-m} = \sum_{n=0}^{\infty} \binom{m-1+n}{m-1} w^n$$
Then
$$\lim_{n \to \infty} \frac{a_n}{a_{n+1}} = \lim_{n \to \infty} \frac{n+1}{n+m} z_0$$
EDIT
@TCL observed that we can simply require that $b_n z_0^n$ goes to zero as $n \to \infty$. Then for a simple pole
$$\lim_{n \to \infty} \frac{a_n}{a_{n+1}} = \lim_{n \to \infty} \frac{\displaystyle b_n z_0^n + \frac{c}{z_0}}{\displaystyle b_{n+1} z_0^n + \frac{c}{z_0^2}}$$
which you can see goes to $z_0$.
a) $|f_{\epsilon}(z) - f(z)| = |\epsilon g(z)| < |f(z)|$ on $\partial D(0,1)$ for $\epsilon$ sufficiently small. So we apply Rouché's theorem and we conclude.
b) We use the generalized argument principle: with the usual hypotheses and $g$ holomorphic $$\frac{1}{2\pi i}\int_{\gamma} \frac{f'(z)}{f(z)}g(z) dz = \sum_{k}g(z_{k}) - g(p_{k})$$ where $z_{k}$ are the zeros of $f$ inside $\gamma$ and $p_{k}$ are the poles.
So in our case we have $$\frac{1}{2\pi i}\int_{\gamma} \frac{f'(z)+\epsilon g'(z)}{f(z) + \epsilon g(z)}z \ dz = z_{\epsilon}$$ where $\gamma = \partial D(0, 1)$
Best Answer
$f$ has no zeros on the boundary of the unit circle $\Bbb D$, therefore $$ \varepsilon_0 = \frac{ \min \{ |f(z)| : |z| = 1 \}}{1 + \max \{ |g(z)| : |z| = 1 \} } \, . $$ is strictly positive. Then for $0 \le \varepsilon \le \varepsilon_0$ and $|z| = 1$, $$ | (f(z) + \varepsilon g(z)) - f(z) | = \varepsilon |g(z)| < |f(z)| $$ and it follows from Rouché's theorem that $f(z) + \varepsilon g(z)$ and $f(z)$ have the same number of zeros in $\Bbb D$. Since $f$ has exactly one (simple) zero, it follows that $f(z) + \varepsilon g(z)$ also has exactly one zero $z_\varepsilon$ in the unit disk.
The solution $z_\varepsilon$ of $f(z) + \varepsilon g(z) = 0$ can be represented as $$ z_\varepsilon = \frac{1}{2 \pi i} \int_{\partial \Bbb D} \frac{z (f'(z) + \varepsilon g'(z))}{f(z) + \varepsilon g(z)} \, dz $$
(Proof: For fixed $\varepsilon$, $f(z) + \varepsilon g(z) = (z - z_\varepsilon)h(z)$ where $h$ is holomorphic and not zero in $\Bbb D$. Then $$ \frac{z (f'(z) + \varepsilon g'(z))}{f(z) + \varepsilon g(z)} = z \frac{h'(z)}{h(z)} + 1 + \frac{z_\varepsilon}{z_\varepsilon - z} $$ and that has exactly one (simple) pole in $\Bbb D$, with residue $z_\varepsilon$.)
Since the integrand is continuous as a function of $(z, \varepsilon) \in \partial \Bbb D \times [0, \varepsilon_0]$ it follows that the integral is a continuous function of $\varepsilon $.
(It is even an analytic function of $\varepsilon $ if we consider complex $\varepsilon $ with $|\varepsilon| < \varepsilon_0 $.)
An alternative proof would be to apply the implicit function theorem to $F(\varepsilon, z) = f(z) + \varepsilon g(z)$, viewed as a function from $\Bbb R \times \Bbb R^2 \to \Bbb R^2$. Writing $f(z) = u(x, y) + i v(x,y)$, the derivative of $F$ with respect to $(x, y)$ at $(0, (0, 0))$ is $$ \begin{pmatrix} u_x(0, 0) & u_y(0, 0) \\ v_x(0, 0) & v_y(0, 0) \end{pmatrix} = \begin{pmatrix} u_x(0, 0) & u_y(0, 0) \\ -u_y(0, 0) & u_x(0, 0) \end{pmatrix} $$ which is invertible because its determinant $$ u_x(0, 0)^2 + u_y(0, 0)^2 = |f'(0)|^2 $$ is not zero.