How to numerically perform analytic continuation? I understand that one can in theory analytically continue a function by repeatedly computing new Taylor series. Suppose for example we have an analytic function $f$ defined on some open set $U$ and compute
$$T_0(z)=\sum_{n=0}^\infty\frac{f^{(n)}(z_0)}{n!}(z-z_0)^n$$
for some $z_0\in U$ close to the boundary of $U$. If this converges on $V$ where $U\cap V$ is non-empty, we can then compute another Taylor series to extend further:
$$T_1(z)=\sum_{n=0}^\infty\frac{T_0^{(n)}(z_1)}{n!}(z-z_1)^n$$
for some $z_1\in V\setminus U$ etc.
However, it is impossible to compute infinitely many terms and higher derivatives quickly become prone to large amounts of cancellation error.
Furthermore, one must repeated drop the degree of the next series expansion, as demonstrated here, in order for the result to be useful. Otherwise, with the same degree at the new point $z_1$, you will end up recovering the original $T_0$ and fail to approximate $f$ further away.

So how can one actually numerically compute the analytic continuation of a function?

In my specific case, I have a set of data points over a subinterval of $\mathbb R$ and I know some basic behavior about the function $f$'s derivatives (all derivatives are positive over the subinterval and to the right, which is the area I want to continue to) and that it has no singularities to the right of the given subinterval.
 A: It's often possible to increase the radius of convergence by applying a suitable conformal mapping to the series. The transform then moves singularities away to a larger distance from the point you are expanding around. For example, the conformal mapping:
$$z = \frac{p w}{p+1-w}\tag{1}$$
Keeps the points $z=0$ and $z=1$ fixed while it moves the point $z = -p$ to infinity, while the point at infinity, which is often a singular point, is moved to $w = p+1$.
Even if only a few terms of the series expansion of a function are known, then we can apply such a conformal transform to evaluate the function outside the radius of convergence. Take e.g. the series for $\log(1+2 z)$:
$$\log(1+2 z) = 2 z - 2 z^2 +\frac{8 z^3}{3}-4 z^4+\frac{32 z^5}{5}-\frac{32 z^6}{3}+\frac{128 z^7}{7} -32 z^8 +\frac{512 z^9}{9} -\frac{512 z^{10}}{5}\cdots$$
This series has a radius of convergence of $\frac{1}{2}$, so it seems that we cannot use this series to evaluate  $\log(1+2 z)$ at $z = 1$. However, we can use this series by applying the conformal mapping (1) for e.g. $p = 1$. This yields the series:
$$w + \frac{w^3}{12} + \frac{w^5}{80} + \frac{w^7}{448} + \frac{w^9}{2304}+\cdots$$
The point $z = 1$ corresponds to $w = 1$, and it's clear that the series converges very fast for $w = 1$, while the original series was divergent at $z = 1$. Note that the computation of this series involves just the substitution of (1) in the series and then a re-expansion in powers of $w$ to order 10.
A: Too longfor a comment.
Using @Saibal Mitra's example
$$\log(1+2z)=2 z-2 z^2+\frac{8 z^3}{3}-4 z^4+\frac{32 z^5}{5}-\frac{32
   z^6}{3}+\frac{128 z^7}{7}+O\left(z^8\right)$$ does not allow to do anything if $z > \frac 12$. For $z=1$, it would give the ridiculous $\frac{444}{35}$ which is $11.5$ times too large.
Using the same coefficients, consider the $[4,3]$ corresponding Padé approximant
$$\log(1+2z)=\frac{2 z+\frac{34 z^2}{7}+\frac{8 z^3}{3}+\frac{4
   z^4}{35} } { 1+\frac{24 z}{7}+\frac{24 z^2}{7}+\frac{32
   z^3}{35}}+O\left(z^8\right)$$  For $z=1$, it gives $\frac{1012}{921}$ which corresponds to a relative error of $0.018$%. For $z=2$, it leads to a relative error of $0.22$%.
