To illustrate the method in it's simplest form lets take the 2D problem without the timing offset as an example. We have the equations
$$\sqrt{(x-x_1)^2+(y-y_1)^2} = d_1$$
$$\sqrt{(x-x_2)^2+(y-y_2)^2} = d_2$$
where $x,y$ are the unknowns. To solve this using Newton's method we rephase it on the form ${\bf F}(x,y) = 0$ for some function ${\bf F}$. A natural choice to take is ${\bf F}(x,y) = (f_1(x,y),f_2(x,y))$ where
$$f_1(x,y) = \sqrt{(x-x_1)^2+(y-y_1)^2}-d_1$$
$$f_2(x,y) = \sqrt{(x-x_2)^2+(y-y_2)^2}-d_2$$
Now Newton's method says that given a starting guess $\vec{r}_0 = (x_0,y_0)$ we can find a new solution by itterating the recurrence
$${\bf J}(\vec{r}_n)(\vec{r}_{n+1}-\vec{r}_n) = -{\bf F}(\vec{r}_n)$$
or equivalently
$$\vec{r}_{n+1} = \vec{r}_n - {\bf J}^{-1}(\vec{r}_n){\bf F}(\vec{r}_n)$$
where ${\bf J}$ is the Jacobian matrix (and ${\bf J}^{-1}$ is the inverse Jacobian)
$${\bf J}(x,y) = \left(\matrix{\frac{\partial f_1}{\partial x} & \frac{\partial f_1}{\partial y}\\ \frac{\partial f_2}{\partial x} & \frac{\partial f_2}{\partial y}}\right) = \left(\matrix{\frac{x-x_1}{\sqrt{(x-x_1)^2+(x-y_1)^2}} & \frac{y-y_1}{\sqrt{(x-x_1)^2+(x-y_1)^2}}\\ \frac{x-x_2}{\sqrt{(x-x_2)^2+(x-y_2)^2}} & \frac{x-y_2}{\sqrt{(x-x_2)^2+(x-y_2)^2}}}\right)$$
and $\vec{r}_n = (x_n,y_n)$. Apposed to Newton's method in 1D we need to solve a linear equation system ($Ax = b$) at every step (unless we can analytically invert the Jacobian) to find the new solution. To demonstrate the method in action, here is a Mathematica implementation:
(* Desired solution *)
xx = 2.4; yy = 4.2;
(* Constants *)
x1 = 0.5; x2 = 2.4; y1 = 0; y2 = -1.3;
d1 = Sqrt[(xx - x1)^2 + (yy - y1)^2];
d2 = Sqrt[(xx - x2)^2 + (yy - y2)^2];
(* F = 0 function *)
f1[x_, y_] = Sqrt[(x - x1)^2 + (y - y1)^2] - d1;
f2[x_, y_] = Sqrt[(x - x2)^2 + (y - y2)^2] - d2;
(* Jacobian matrix *)
J[x_, y_] = {{D[f1[x, y], x], D[f1[x, y], y]}, {D[f2[x, y], x], D[f2[x, y], y]}};
(* Random first guess *)
rn = {-3.7, 8.9};
Do[
(* Itterate *)
rn += -Inverse[J[rn[[1]], rn[[2]]]].{f1[rn[[1]], rn[[2]]], f2[rn[[1]], rn[[2]]]};
Print["Step ",i," Solution = ", rn];
, {i, 1, 10}];
Here is a test run. After only $6$ steps we have found the true solution to high accuracy:
Step 1 Solution = $\{9.4199,9.30668\}$
Step 2 Solution = $\{-0.0657397,6.9274\}$
Step 3 Solution = $\{3.95296,4.90711\}$
Step 4 Solution = $\{2.24599,4.40806\}$
Step 5 Solution = $\{2.4086,4.20223\}$
Step 6 Solution = $\{2.4,4.20001\}$
Step 7 Solution = $\{2.4,4.2\}$
I did not add a convergence criterion above, usually one can do this by defining convergence when $|{\bf F}(\vec{r}_n)| < \epsilon$ where $\epsilon$ is some small pre-defined constant.
To generalize this procedure to your problem you need to add two more variables ($z$ and ${\rm d}T$) and add two more functions to ${\bf F}$ and then calculate the now $4\times 4$ Jacobian matrix.
Best Answer
Let be $e_n=\left\vert x_{n+1}-x_n\right\vert$. You can estimate the value of $m$ considering $$e_{n+1}=\frac{m-1}{m}e_n.$$ Just make two normal Newton steps and use $x_0,x_1,x_2$ to evaluate the two succesive errors. The first steps of any iterative method may not be very good, specially in multiple roots, so you can discard $x_0,x_1$ and take two additional steps before approximate $m$.
Deepening a little more: if a function $f$ has a multiple root $r$ with multiplicity $m$, then $f$ is of the form $$f(x)=(x-r)^mg(x)$$ with necessarly $g(r)\neq 0$ and $f^\prime(x)=f^{\prime\prime}(x)=\cdots =f^{(m-1)}(x)=0$. Let be $u(x)=f(x)/f^\prime (x)$. The function $u(x)$ has only one root (easy to proof) and you can use the iteration $$x_{n+1}=x_n-\frac{u(x)}{u^\prime (x)}$$ and it converges cuadratically.
Also you can dynamically estimate $m$ with $$m=\max\left\{1, INT\left(\frac{x_n-x_{n-1}}{u(x_n)-u(x_{n-1})}\right)\right\}$$ as William Kahan suggest.