How to solve this system of hyperbola equations? I have a system of 2 equations each describing the branch of a hyperbola. The below equations represent hyperbolae with foci $P_0$ and $C_1$ (or $C_2)$ and transverse axis length $r_1$ (or $r_2)$.  I'm trying to find the unknowns $x$ and $y$:
$$\sqrt{(x-x_1)^2+(y-y_1)^2}-\sqrt{(x-x_0)^2+(y-y_0)^2}=r_1. \tag{1}$$
$$\sqrt{(x-x_2)^2+(y-y_2)^2}-\sqrt{(x-x_0)^2+(y-y_0)^2}=r_2. \tag{2}$$
Where:
$P_0 = (x_0, y_0) = (0.87,-0.5)$
$C_1 = (x_1, y_1, r_1) = (0, 1, 0.13)$
$C_2 = (x_2, y_2, r_2) = (-0.87, -0.5, 0.49)$
These hyperbolae are the result of an algorithm for trilateration I'm trying to compute, known as Time Difference Of Arrival (TDOA). I just have no idea where to go from here!
EDIT: For some additional context, here's a graph of the problem:

 A: 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$.
A: The system of equations corresponds to a special case CCP of Apollonius' problem. Namely, $C=(x,y)$ is the center of a circle $\overline{C}$ with radius $r=\sqrt{(x-x_0)^2+(y-y_0)^2}$, which is the outer tangent to a circle $\overline{C_1}$ with the center $C_1$ and radius $r_1$, to a circle $\overline{C_2}$ with the center $C_2$ and radius $r_2$, and passing through the point $P_0=(x_0,y_0)$. Let $a=|P_0C_1|=|P_0C_2|=|C_1C_2|$. Let $i=1$ or $i=2$. 
If $r_i>a$ then $P_0$ belongs to the inner region of the circle $\overline{C_i}$, so the circle $\overline{C}$ passing through $P_0$
cannot be outer tangent to the circle $\overline{C}$. 
If $r_i=a$ then the circle $\overline{C_i}$ passes through the point $P_0$. Let $\ell_i$ be the tangent line to $\overline{C_i}$ at $P_0$. Then the circle $\overline{C}$ is tangent to $\overline{C_i}$ at $P_0$ and lies at the other halfplane with respect to $\ell_i$ than $\overline{C_i}$. Thus $C$ lies on the ray $r_0$, emanating from $P_0$ and extending the segment $C_iP_0$. Let $i’=3-i$. The circle $\overline{C}$ can touch the circle $\overline{C_{i’}}$ ony if $r_{i’}>d$, where $d=a\sin\tfrac{\pi}{6}=\tfrac{a}{2}$ is the distance from $C_{i’}$ to $\ell_i$. On the other hand, it is easy to see that if $r_{i’}>d$ then $\overline{C}$ touches $\overline{C_{i’}}$ iff $(r+r_{i’})^2=(r+d)^2+\left(a\cos\tfrac{\pi}{6}\right)^2$, that is iff $r=(a^2-r_{i’}^2)(2r_{i’}-a)$.
So farther we assume that $r_i<a$. In order to find $(x,y)$ we shall follow the construction proposed in Wikipedia. Make an inversion centered at $P_0$ and the radius $a$ . For $i=1$ or $i=2$ let $\overline{C_i}’$ be the circle in which the inversion transforms the circle $\overline{C_i}$. The inversion transforms $\overline{C}$ into a line $\overline{C}’$ which is a common tangent to $\overline{C_1}’$  and $\overline{C_2}’$. Let $n’$ be a unit (that is of length $1$) normal vector for $\overline{C}’$ and $z’\in \overline{C}’$ be a point such that $z’=P_0+\lambda n’$ for some $\lambda>0$. That is $\lambda$ is the distance from the point $P_0$ to the line $\overline{C}’$. Let $z=P_0+\tfrac {a^2}{\lambda} n’$ be the inversion image of $z’$. Then $P_0$ and $z$ are diametral points of the circle $\overline{C}$, so its center $(x,y)$ is $\tfrac 12(P_0+z)=P_0+\tfrac {a^2}{2\lambda} n’$. So in order to find $(x,y)$ it suffices to find $n’$ and $\lambda$.
Let $z_i=C_i-P_0=(x_i,y_i)-(x_0,y_0)$. Clearly, $|P_0C_i|=a$. Let $s_i=P_0+z_i\tfrac{a+r_i}{a}$ and $t_i=P_0+z_i\tfrac{a-r_i}{a}$ be the intersection points of the circle $\overline{C_i}$ with the line spanned by the points $P_0$ and $C_i$.  Let $s’_i=P_0+z_i\tfrac{a}{a+r_i}$ and $t’_i=P_0+z_i\tfrac{a}{a-r_i}$ be the inversion images of the points $s_i$ and $t_i$, respectively. Then $t’_i$ and $s’_i$ are diametral points of the circle $\overline{C_i}’$. So it has a center $C_i^*=\tfrac 12(s_i’+t_i’)=P_0+z_i\tfrac{a^2}{a^2-r_i^2}$ and a radius $r_i^*=\tfrac 12|s_i’-t_i’|=\tfrac {ar_i^2}{a^2-r_i^2}$. Since the common tangent $\overline{C}’$ is perpendicular to $n’$, contains circles $\overline{C_1}’$, $\overline{C_2}’$, and the point $P_0$ is one of its halfplanes,  considering orthogonal projections of $C_i^*$ onto $n’$, we obtain $\lambda=(C_i^*-P_0, n’)+r_i^*$.  Put $C^*= C_2^*-C_1^*$ and $r^*=r_2^*-r_1^*$. Then $(C^*, n’)+r^*=0$.
Suppose to the contrary that $C^*=0$. Then $r^*=0$, that is $C_2^*=C_1^*$ and $r_2^*=r_1^*$. Then the circles $\overline{C_1}’$ and $\overline{C_2}’$ coincide, and the circles $\overline{C_1}$ and $\overline{C_2}’$ coincide too, which is imposible becuse their centers are distinct. So $C^*\ne 0$. 
Let $0\le \varphi<2\pi$ be the angle between the vectors $C^*$ and $n’$. Since $(C^*, n’)=|C^*||n’|\cos\varphi$ and $n’$ is a unit vector, we see that there are at most two possible values for $\varphi$, and so for $n’$. 
Assume that there are two solutions $n_1$ and $n_2$ such that $\varphi_j$ is the angle between the vectors $C^*$ and $n’_j$ for each $j$ then we have $\varphi_1+\varphi_2=2\pi$. This follows that $n’_1+n’_2=\mu C^*$ for some real $\mu$. 
Since we already know $C_i^*$ and $r_i^*$,  $(C^*, n’)+r^*=0$ is a linear equation for conponents $n_x’$ and $n_y’$ of the vector $n’$. Since $n’$ is a unit vector, we also have  $(n_x’)^2+(n_y’)^2=1$. So we can find both possible values of $n’$ solving a quadratic equation. Then we calculate $\lambda=(C_i^*-P_0, n’)+r_i^*$ and $(x,y)=P_0+\tfrac {a^2}{2\lambda} n’$.
