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:

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$.