As per a comment of Hagen von Eitzen: the complement of two disjoint disks is not simply connected. But that region can be mapped to an annulus.
Assume two circles, $C_1$ and $C_2$ with radius $r_1$ and $r_2$, respectively, such that $r_1 > r_2$. Assume yet that $C_1$ is centered at the origin and $C_2$ centered over the Real axis at a point $h$. Lets transform $C_1$ and $C_2$ such that they form an annulus centered at the origin.
The algorithm is as follows:
Find parameters $a,b,c,d$ of the Möbius transformation (can be by explicit determinant formula) that makes the bigger circle into a line over the Imaginary axis ($\mathbf{z} \to \mathbf{w}$).
$$
\mathbf{z} = [i\,r_1,\, -i\,r_1,\, r_1] \\
\mathbf{w} = [-i,\, i,\, \infty] \\
$$
This transformation is
$$
f(z) = \frac{a\,z + b}{c\,z + d}
$$
Move both the line and new circle by $-p$ such that the circle of radius $r_w$ is centered at the origin. Variables $p$ and $r_w$ satisfy:
$$
|w_1|^2 - 2\,\Re(w_1 \,\bar p) + |p|^2 - r_w^2 = 0 \\
|w_2|^2 - 2\,\Re(w_2 \,\bar p) + |p|^2 - r_w^2 = 0 \\
$$
Where $w_1$ and $w_2$ are points over the circle.
This transformation is
$$
g(z) = z - p
$$
Now find $x_1$ and $x_2$ (zeros of the following equation) that transforms the line into a circle around the origin while also preserving the center of the other circle around the origin.
$$
x^2 + 2\,p\,x + r_w^2 = 0
$$
This transformation is
$$
h(z) = \frac{z - x1}{z - x2}
$$
Now compose all these transformations:
$$
T(z) = h \circ g \circ f (z)
$$
$$
T(z) = \frac{b + a\,z - (p + x_1) (d + c\,z)}{b + a\,z - (p + x_2) (d + c\,z)}
$$
A Julia code
The following is a Julia code I created that implements the above algorithm.
using LinearAlgebra
using NLsolve
using Polynomials
using Plots
function mobius(z, w)
a = det([z[1]*w[1] w[1] 1;
z[2]*w[2] w[2] 1;
z[3]*w[3] w[3] 1])
b = det([z[1]*w[1] z[1] w[1];
z[2]*w[2] z[2] w[2];
z[3]*w[3] z[3] w[3]])
c = det([z[1] w[1] 1;
z[2] w[2] 1;
z[3] w[3] 1])
d = det([z[1]*w[1] z[1] 1;
z[2]*w[2] z[2] 1;
z[3]*w[3] z[3] 1])
return a, b, c, d
end
function find_trans(r1, r2, h, n=200)
infty = 1e10 # "infinity"
θ = range(0, 2π, length=n)
z1 = @. r1*(cos(θ) + 1im*sin(θ))
z2 = @. r2*(cos(θ) + 1im*sin(θ)) + h
# transforms the bigger circle into a line on the imaginary axis
a, b, c, d = mobius([-1im*r1, 1im*r1, r1], [1im, -1im, infty])
f(z) = (a*z + b)/(c*z + d)
# center transformed circle on the origin
w1 = f(z2[1])
w2 = f(z2[cld(n,2)])
function f!(F, x)
F[1] = abs(w1)^2 - 2*real(w1*conj(x[1])) + abs(x[1])^2 - x[2]^2
F[2] = abs(w2)^2 - 2*real(w2*conj(x[1])) + abs(x[1])^2 - x[2]^2
end
# sol.zero: center point of the transformed circle and its radius
sol1 = nlsolve(f!, [0.0, 2r1])
p = sol1.zero[1]
r = sol1.zero[2]
# transform the line into a unit circle centered at the origin,
# while also preserving the other at the origin
x1, x2 = roots(Poly([r^2, 2p, 1]))
# Möbius transformation that sends x1 -> 0 and x2 -> oo,
# where x1 and x2 are symmetric points to both the line and circle
T(z) = (b + a*z - (p + x1)*(d + c*z))/(b + a*z - (p + x2)*(d + c*z))
return a, b, c, d, p, r, x1, x2, T
end
begin
n = 500
θ = range(0, 2π, length=n)
r1 = 2
r2 = 1
h = 2r1
z1 = @. r1*(cos(θ) + 1im*sin(θ))
z2 = @. r2*(cos(θ) + 1im*sin(θ)) + h
a, b, c, d, p, r, x1, x2, T = find_trans(r1, r2, h)
w1 = T.(z1)
w2 = T.(z2)
plot(aspect_ratio=1)#, xlim=(-l,l), ylim=(-l,l))
plot!(z1, color=:blue, label="z1")
scatter!(z1[1:1], color=:blue, label="")
plot!(z2, color=:red, label="z2")
scatter!(z2[1:1], color=:red, label="")
plot!(w1, color=:blue, label="w1", linestyle=:dash)
scatter!(w1[1:1], color=:blue, label="")
plot!(w2, color=:red, label="w2", linestyle=:dash)
scatter!(w2[1:1], color=:red, label="")
end
The output of the above code is the following figure.
Now the potential problem can be solved in the annulus. See, for example, Conformal Mapping and Bipolar Coordinate for Eccentric Laplace Problems.
Edit
For the sake of completness, the inverse transformation is
$$
T^{-1}(z) = \frac{b - d\,(p + x_1) - b\,z + d\,(p + x_2)\,z}{a\,(z - 1) + c\,(p + x_1 - (p + x_2)\,z)}
$$
