Reducing this non-linear system I'm trying to solve the following non-linear system:
$(A_x - P_x)^2 + (A_y - P_y)^2 + (A_z - P_z)^2 = (v(t_a - t_0))^2$ $\{1\}$
$(B_x - P_x)^2 + (B_y - P_y)^2 + (B_z - P_z)^2 = (v(t_b - t_0))^2$ $\{2\}$
$(C_x - P_x)^2 + (C_y - P_y)^2 + (C_z - P_z)^2 = (v(t_c - t_0))^2$ $\{3\}$
$(D_x - P_x)^2 + (D_y - P_y)^2 + (D_z - P_z)^2 = (v(t_d - t_0))^2$ $\{4\}$
Knowns:
$A_x, A_y, B_x, B_y, C_x, C_y, D_x, D_y, v, t_a, t_b, t_c, t_d$
Unknowns: $P_x, P_y, t_0$
The end goal is to determine $P_x, P_y$

I know that a slight variation of this system can be solved in the following manner.
$(A_x - P_x)^2 + (A_y - P_y)^2 = (v(t_a - t_0))^2$ $\{1\}$
$(B_x - P_x)^2 + (B_y - P_y)^2 = (v(t_b - t_0))^2$ $\{2\}$
$(C_x - P_x)^2 + (C_y - P_y)^2 = (v(t_c - t_0))^2$ $\{3\}$
Now, open the brackets, and form two new equations, {2'} = {3}-{2}, and {3'} = {1} - {3}:
$C_x^2-B_x^2 -2*(C_x-B_x)*P_x + C_y^2-B_y^2 -2*(C_y-B_y)*P_y = v^2*(t_c^2-t_b^2 -2*(t_c-t_b)*t_0)$
$A_x^2-C_x^2 -2*(A_x-C_x)*P_x + A_y^2-C_y^2 -2*(A_y-C_y)*P_y = v^2*(t_a^2-t_c^2 -2*(t_a-t_c)*t_0)$
These equations have the form:
$HP_x + JP_y = Kt_0 + L$
$TP_x + UP_y = Vt_0 + W$
Calculating $T(2') - H(3')$ and $U(2') - J(3')$ gives:
$P_x = \frac{T(Kt_0 + L) - H(Vt_o + W)}{TJ-HU}$, $P_y = \cdots$
Substitute these into {1} to obtain a quadratic equation in $t_0$. We then simple solve for the root of this quadratic, and plug into the formulae for $P_x$, $P_y$.

I'd like to follow a similar process for this system. I'll set aside {1} (in the first system), and form the linear system {2'} = {2}-{1}, {3'} = {3} - {1}, and {4'} = {4} - {1}. With these, I should be able to do the same as what I'd done above, obtaining an expression for $P_x, P_y, P_z$ in terms of $t_0$.
I then substitute these into {1} to obtain a quadratic in terms of $t_0$. I solve for $t_0$, then plug into the formulae to find $P_x, P_y, P_z$
My problem, and what I've tried:
I'm not sure what precisely this would look like. I've spent some time attempting to perform the manipulations by hand to obtain the expressions for $P_x, P_y, P_z$ and the quadratic in $t_0$, however I've failed. I've gotten several different answers, none of which look right to me.
Moreover, I've tried running this through a computer algebra system, and it choked up (unfortunately, I do not have access to something "more popular" like Mathematica). I tried running it on my school's CAS software, and it too failed. Even my maths professor has come up with multiple answers.
What should $P_x, P_y, P_z$ and the quadratic look like?

Note: I'm aware that I'll end up with two roots from the quadratic. In my case, this is ok, I want both.

 A: You should not compute anything directly with components. Instead, you should express the relation in terms of matrices. There are enough structure among your variables and finding the quadratic equation is more or less equivalent to inverting some $4 \times 4$ matrix.
Instead of points in $3$-d, think in terms of events.
For each event, say point $A = (A_x,A_y,A_z)$ at time $t_A$, we will associate a $4\times 1$ column vector $(A_x,A_y,A_z,v t_a)^T$ to it. Let

*

*$u_1, u_2, u_3, u_4$ be the column vectors corresponding to $A, B, C, D$.

*$p$ be the column vector $(P_x,P_y,P_z,vt)^T$.

*$\eta$ be the $4 \times 4$ diagonal matrix with entries $(1,1,1,-1)$ along
its diagonal.

In terms of these, you equations becomes
$$(p - u_i)^T \eta (p - u_i) = 0
\iff 2 u_i^T \eta p = p^T\eta p + u_i^T\eta u_i
\quad\text{ for } 1 \le i \le 4
$$
Let $s = p^T\eta p$ and

*

*$\Omega$ be the $4 \times 4$ matrix whose $i^{th}$ row is $2u_i^T$.

*$U$ be the $4 \times 1$ column vector whose $i^{th}$ entry is $u_i^T\eta u_i$.

*$\theta$ be the $4 \times 1$ column vector with all entries $1$.

In terms of these, your set of equations become
$$\Omega \eta p = s \theta + U
\implies p = \eta \Omega^{-1}(s \theta + U)\tag{*1}
$$
Substitute this back into definition of $s$, one obtain a quadratic equation in $s$.
$$s = p^T \eta p = (s\theta + U)^T\Omega^{-T}\eta \Omega^{-1}(s \theta + U)$$
Solving this give you two possible values of $s$. Plug that into $(*1)$ will give you the $p$ and hence $(P_x,P_y,P_z,t)$.
As one can see, if one doesn't insists writing down everything in terms of components, setting up the column vectors $u_i, U$ and matrices $\Omega$
is pretty straight forward. The only thing that is slightly more complicated is taking the inverse of the $4\times 4$ matrix $\Omega$.
