As you said, the plane that you're looking for has the equation
$ n \cdot (x - x_0) = 0 $
The plane has to satisfy that the distance of the centers from the plane are $r_1, r_2, r_3$.
So if $x_1 , x_2, x_3$ are the position vectors of the three centers then we want to solve
$ | n \cdot (x_1 - x_0) | = r_1 $
$ | n \cdot (x_2 - x_0) | = r_2 $
$ | n \cdot (x_3 - x_0) | = r_3 $
where $n$ also satisfies, $n \cdot n = 1 $
To solve these four equations, define
$ s_1 = \operatorname{sgn}(n \cdot (x_1 - x_0) ) $
$ s_2 = \operatorname{sgn}( n \cdot (x_2 - x_0) ) $
$ s_3 = \operatorname{sgn}(n \cdot (x_3 - x_0 ) ) $
where $\operatorname{sgn}$ is the sign function returning values of $+1$ or $-1$.
Then, assuming $r_1 , r_2 , r_3 \gt 0 $, the equations become
$ n \cdot (x_1 - x_0) = s_1 r_1 $
$ n \cdot (x_2 - x_0) = s_2 r_2 $
$ n \cdot (x_3 - x_0) = s_3 r_3 $
Taking the difference between the first and second equations, and then the difference between the first and third equations, leads to
$ n \cdot (x_1 - x_2) = s_1 r_1 - s_2 r_2 $
$ n \cdot (x_1 - x_3) = s_1 r_1 - s_3 r_3 $
This is a linear system of $2$ equations in $3$ unknowns (which are the components of $n$). Therefore, the solution will take the form
$ n = n_0 + t n_1 $
where $n_0$ and $n_1$ have been obtained by the Gauss-Jordan elimination process, and $t \in \mathbb{R}$ is yet to be determined.
To determine $t$, we just impose $n \cdot n= 1 $, and this translates into,
$ (n_1 \cdot n_1) t^2 + 2 t (n_1 \cdot n_0) + (n_0 \cdot n_0) - 1 = 0 $
This will give $0, 1,$ or $2$ solutions. Once we have $n$ we can obtain the constant $(n \cdot x_0)$ from the first of the original equations, namely,
$ d = n \cdot x_0 = n \cdot x_1 - s_1 r_1 $
so that the plane is completely determined.
Note that if we reverse all the signs of the linear system, then the obtained normal vector will be just the negative of the one obtained with the original signs. So if we obtain $ n $ with $s_1 = 1 , s_2 = -1, s_3 = -1 $, then we will obtain $ - n $ with $s_1 = -1, s_2 = 1, s_3 = 1 $. Therefore, we only need to take the following
combinations
$ ( +, + , + ) , (+, +, - ) , (+, - , + ), ( -, + , + ) $
And this will generate two planes for each of the signs triple, making a total of $8$ possible planes.
Best Answer
Rewrite the system of equations in their original form,
$$(x-x_0)^2+(y-y_0)^2 =r^2\tag{1}$$ $$(x-x_1)^2+(x-y_1)^2 = (r_1+r)^2\tag{2}$$ $$(x-x_2)^2+(x-y_2)^2 = (r_2+r)^2\tag{3}$$
where $r$ is the radius of the circle with the center $(x,y)$ to be solved.
As seen below, $x$ and $y$ depend on $r$ linearly. It is more convenient to maintain symmetry by solving for $r$ first and then for $(x,y)$.
From (2)-(1) and (3)-(1), we get a pair of linear equations,
$$a_1x+b_1y=c_1+2r_1r\tag{4}$$ $$a_2x+b_2y=c_2+2r_2r\tag{5}$$
where the coefficients are all known and given by,
$$a_1=2(x_0-x_1), \>\>\>\> b_1=2(y_0-y_1),\>\>\>\>\> c_1=r_1^2+x_0^2+y_0^2-x_1^2-y_1^2$$ $$a_2=2(x_0-x_2), \>\>\>\> b_2=2(y_0-y_2),\>\>\>\> c_2=r_2^2+x_0^2+y_0^2-x_2^2-y_2^2$$
Solve the linear equations (4) and (5) for $x$ and $y$ in terms of $r$,
$$x=d_1r+e_1,\>\>\>\>\> y = d_2r+e_2\tag{6}$$
where the known coefficients are,
$$d_1 = \frac{2r_1b_2-2r_2b_1}{a_1b_2-a_2b_1},\>\>\>\>\>\>e_1 = \frac{c_1b_2-c_2b_1}{a_1b_2-a_2b_1} $$
$$d_2 = \frac{2r_1a_2-2r_2a_1}{a_2b_1-a_1b_2},\>\>\>\>\>\>e_2 = \frac{c_1a_2-c_2a_1}{a_2b_1-a_1b_2} $$
Plug (6) back into (1) to get a quadratic equation in $r$,
$$ar^2+br+c=0\tag{7}$$
with the coefficients given by,
$$a=d_1^2+d_2^2-1,\>\>\> b= 2(e_1-x_0)d_1+2(e_2-y_0)d_2,\>\>\>c=(e_1-x_0)^2+(e_2-y_0)^2$$
The equation (7) for $r$ can be solved with the standard quadratic formula,
$$r = \frac{-b\pm \sqrt{b^2-4ac}}{2a}$$
Afterwards, plug $r$ into (6) to obtain the solutions for $x$ and $y$.