(Expanding my comment on James' answer)
James already provided a nice picture, which I'll build upon here:
At the moment, I don't know of an analytical formula, but the algorithm is this:
- Given: the Laurent-series in the yellow zone, delimited by $r<|z-z_0|<R$
- Sought: the one in the next-outer blue zone, delimited by $R<|z-z_0|<{\color{blue}{R_b}}$ (it's straightforward to modify all this for switching to a smaller annulus)
- Shift the centre to a point $z_1$ which lies not in direction of a singularity that causes $R$ to be $R$ (which can be verified by checking the new $R_1$)
- Use the $z_1$-series to obtain values in a small section of the blue annulus around $\color{lime}{z_3}$, put them into the formula for a Laurent-series expansion (actually, Taylor, since we avoided the singularity - to do expand around the singularity instead for fewer steps?) in that small section of blue, with convergence radius $\color{lime}{R_2}$
- move along the yellow-blue interface towards the convergence boundary at $\color{orange}{z_3}$, Taylor-expand there
- go on until a full circle is available (when approaching singularities on the interface, go round them via the blue zone)
- use the Laurent-expansion with that bumped-circle around $z_0$
Ok, now let's try to get this more analytical:
The series representation around $z_l$ is
\begin{align}
f_l(z) &= \sum_{k=-\infty}^\infty c_k^{(l)} (z-z_l)^k = f(z)\Big|_{r_l<|z-z_l|<R_l},
\\ &\quad\text{and in its convergence zone we have}
\\ c_k^{(l)} &= \frac1{2\pi i}\oint\limits_{r_l<|z-z_l|<R_l}\frac{f(z)\,dz}{(z - z_l)^{k+1}} \quad \Big| \quad z = z_l + re^{i\phi},\ dz = ir e^{i\phi}\,d\phi,\ r\in(r_l,R_l)
\\ &= \frac1{2\pi r^k}\int_0^{2\pi}f(z_l+re^{i\phi})\cdot e^{-ik\phi}\,d\phi. \tag{c}\label{c}
\end{align}
(Note how $\eqref{c}$ is a Fourier transform along a circle, see this question for more on that)
Now, let's use a slightly different centre $z_m = z_l - d_m$ and choose $r\in(r_l+|d_m|,R_l-|d_m|)$ (which implicitly requires $|d_m|<\frac{R_l-r_l}2$ to make sense) so we remain in the yellow annulus:
\begin{align}
c_k^{(l)} &= \frac1{2\pi r^k}\int_0^{2\pi}\Big[\sum_p c_p^{(m)}\underbrace{(re^{i\phi} + d_m)^p}_{=\sum\limits_{n=0}^p \binom pn d_m^{p-n}r^ne^{in\phi}}\Big]\cdot e^{-ik\phi}\,d\phi \quad\Bigg|\quad \int_0^{2\pi}e^{i(n-k)\phi}\,d\phi = 2\pi\delta_{nk}
\\ &= \sum_{p=k}^\infty \binom pk c_p^{(m)}d_m^{p-k}.
\end{align}
Let's take $d_m$ to be infinitesimal, so we have
$$c_k^{(l)} \dot= c_k^{(m)} + (k+1)c_{k+1}^{(m)}d_m$$
where $\dot=$ denotes equal up to $\mathcal O(d_m^2)$.
Ok, so the next step is obtaining the new convergence radii to make sure we didn't accidentally move towards a singularity:
\begin{align}
\frac1{R^{(l)}} &= \limsup_{k\to\infty}|c_k^{(l)}|^{\frac1k}
\\ &= \limsup_{k\to\infty}\underbrace{\Big|c_k^{(m)} + (k+1)c_{k+1}^{(m)}d_m\Big|^{\frac1k}}_{\dot=\Big(|c_k^{(m)}|^2 + 2(k+1)\Re\big[(c_k^{(m)})^*c_{k+1}^{(m)}d_m\big] \Big)^{\frac1{2k}}}
\\ &\dot= \limsup_{k\to\infty} |c_k^{(m)}|^{\frac1k}\Bigg(1+\tfrac{k+1}{k}\underbrace{\frac{\Re\big[(c_k^{(m)})^*c_{k+1}^{(m)}d_m\big]}{|c_k^{(m)}|^2}}_{=\Re\frac{c_{k+1}^{(m)}d_m}{c_k^{(m)}}}\Bigg)
\\ &\le \frac1{R^{(m)}}\Big(1+\Re\Big[d_m \limsup_{k\to\infty}\tfrac{k+1}k\tfrac{c_{k+1}^{(m)}}{c_k^{(m)}}\Big]\Big)
\\\Rightarrow R^{(l)} & \ge R^{(m)}\underbrace{\Big(1+\Re\Big[d_m \limsup_{k\to\infty}\tfrac{k+1}k\tfrac{c_{k+1}^{(m)}}{c_k^{(m)}}\Big]\Big)^{-1}}_{\dot=1-\Re\Big[d_m \limsup_{k\to\infty}\tfrac{k+1}k\tfrac{c_{k+1}^{(m)}}{c_k^{(m)}}\Big]}
\end{align}
I'm sure that can be improved somehow...
To be continued (next step: calculate new convergence radius to make sure one doesn't approach the interface-singularities)
