What type of series converges around two points (instead of one like in Taylor series)? A converging Taylor series is known to converge in a disk around point of expansion $z_0$. Its partial sum of order $n$ approximates the function so that its error is
$$f_n(x)-f_\infty(x)=\mathcal O((z-z_0)^{n+1}).$$
Also, the function $f_\infty$ to which the series converges is always analytical at $z=z_0$.
I wonder whether it's possible to construct a series which would have two points $z_1$ and $z_2$, for which the error in partial sums would be like
$$g_n(x)-g_\infty(x)=\mathcal O\left((z-z_1)^{n+1}(z-z_2)^{n+1}\right)$$
or something similar, and its convergent $g_\infty$ would be analytical at both $z_1$ and $z_2$. Is there such a thing?
 A: TLDR: yes, you can expand in $(x^2-1)^n$ polynomials. See the procedure around equation $(8)$.

Expansion in orthogonal polynomials like e.g. Legendre polynomials reduces to Taylor expansion when we let the length of domain of expansion $a\to0$. This works for the families of orthogonal polynomials which have dot product defined as a single integral over a line segment.
But we can also generate a family of polynomials which would be orthogonal with respect to another dot product, e.g.
$$\langle f,g\rangle=\int\limits_{-p-a/2}^{-p+a/2} f(x)g(x)\,\mathrm dx+\int\limits_{p-a/2}^{p+a/2} f(x)g(x)\,\mathrm dx.
\tag1$$
Here $p$ is the center of one of the half-domains of convergence, with the other on the opposite side with respect to $x=0$, and $a$ is the length of these half-domains.
We can then start with a set of monomials $\{x^n\}_{n=0}^{\infty}$, and orthogonalize them using e.g. Gram-Schmidt process with respect to dot product $(1)$. E.g. first 3 polynomials would be
$$\begin{align}
f_0&=\frac1{\sqrt{2a}},\\
f_1&=x\sqrt{\frac6{a^3+12ap^2}},\\
f_2&=\frac{\sqrt5(a^2+12(p^2-x^2))}{2\sqrt{2a^5+120a^3p^2}}.
\end{align}
\tag2$$
Then we can expand in this set of polynomials the function we're interested in, and then send $a\to0$. Thus we'll obtain Taylor-like expansion of this function, which will converge around the points $x=\pm p$, similarly to Taylor expansion, which converges around its single point of expansion.

WLOG, let $p=1^\dagger$. Then, if you do the expansion and take the limit as described above, you'll find that you simply get an expansion in the following polynomials:
$$r_n(x)=(x^2-1)^{\lfloor n/2\rfloor}x^{\operatorname{mod}(n,2)},
\tag3$$
where $\operatorname{mod}(n,2)$ evaluates to $1$ whenever $n$ is odd and to $0$ otherwise. Moreover, if you recognize that any function can be represented by a combination of two even functions like
$$f(x)\equiv
\underbrace{\frac{f(x)+f(-x)}2}_\text{even}+
x\cdot\underbrace{\frac{f(x)-f(-x)}{2x}}_\text{even},
\tag4$$
it becomes especially easy to do the expansion of each even component in the polynomials of
$$s_n(x)=(x^2-1)^n.
\tag5$$
$n$th coefficient $b_n$ of expansion of an even function $g$ can be computed as
$$\begin{align}
b_n&=\lim\limits_{x\to1}\frac{g(x)-\sum_{k=0}^{n-1}b_k s_k(x)}{s_n(x)}=\\
&=\lim\limits_{x\to1}\frac{g(x)-\sum_{k=0}^{n-1}b_k (x^2-1)^k}{(x^2-1)^n}.
\end{align}
\tag6$$
Equation $(6)$ is very similar to a way of calculating Taylor coefficients for a function $g$ expansion at $x=1$:
$$t_n=\lim\limits_{x\to1}\frac{g(x)-\sum_{k=0}^{n-1}t_k (x-1)^k}{(x-1)^n}.
\tag7$$
Both limits, $(6)$ and $(7)$, being $0/0$ indeterminate forms, allow straightforward calculation using l'Hôpital's rule $n$ times. Doing this with $(6)$ yields:
$$\begin{align}
b_n&=\frac{\frac{\mathrm d^n}{\mathrm d x^n}\left.\left(g(x)-\sum_{k=0}^{n-1}b_k s_k(x)\right)\right|_{x=1}}{s_n^{(n)}(1)}=\\
&=\frac{g^{(n)}(1)-\sum_{k=0}^{n-1}b_k s_k^{(n)}(1)}{s_n^{(n)}(1)}.
\end{align}
\tag8$$
And we get the expansion for the even function $g$:
$$g(x)=\sum\limits_{n=0}^\infty b_n s_n(x)=\sum\limits_{n=0}^\infty b_n (x^2-1)^n.
\tag9$$
After combining the two even parts as per $(4)$ we can get the final expansion of our general (not necessarily even or odd) function $f$.

Below is an example of expansion of $x\mapsto\sin(x)+\cos(x)$ function around $x=\pm3$. Top image here is the function and its partial-sum approximations, bottom is absolute difference in log scale.

$^\dagger$ We can always rescale function domain to accommodate the fixed points of expansion, similarly to how we do with orthogonal polynomials
