Taylor Series with WLS

I was doing an estimate position from unknown coordinate transmitter using a Taylor series that was rolled out with Weighted Least Squares, but I always found an unwanted results, which was due to the too large 'dox ', ' doy ', and ' doz ' values. How can I find a good approximation position with this method?

' dox ', ' doy ', and ' doz ' are used to update the position

True coordinate Transmitter (xe= 6.5; ye= 0; ze= 2.5;)

clear all; clc;
tic
%receiver
x1= 1; y1= 4.5; z1= 7;
x2= 1; y2= 1; z2= 2.5;
x3= 5.5; y3= 4.5; z3= 3;
x4= 1; y4= 7; z4= 4.5;
%transmitter
xe= 6.5; ye= 0; ze= 2.5;
x= sym('x'); y= sym('y'); z= sym('z');
dox= sym('dox'); doy= sym('doy'); doz= sym('doz');
%distance
r1= sqrt((x-x1)^2+(y-y1)^2+(z-z1)^2);
r2= sqrt((x-x2)^2+(y-y2)^2+(z-z2)^2);
r3= sqrt((x-x3)^2+(y-y3)^2+(z-z3)^2);
r4= sqrt((x-x4)^2+(y-y4)^2+(z-z4)^2);
c=10;
t1= 0.28;
t2= 0.38;
t3= 0.06;
d21= t1*c;
d31= t2*c;
d41= t3*c;
eps21= 0.543;
eps31= 0.4783;
eps41= 0.33802;
f1= d21-eps21;
f2= d31-eps31;
f3= d41-eps41;
%expansion
k= 50;
for i= 1:k
    d1(i)= sqrt((xe(i)-x1)^2+(ye(i)-y1)^2+(ze(i)-z1)^2);
    d2(i)= sqrt((xe(i)-x2)^2+(ye(i)-y2)^2+(ze(i)-z2)^2);
    d3(i)= sqrt((xe(i)-x3)^2+(ye(i)-y3)^2+(ze(i)-z3)^2);
    d4(i)= sqrt((xe(i)-x4)^2+(ye(i)-y4)^2+(ze(i)-z4)^2);
    a11(i)= ((x1-xe(i))/d1(i))-((x2-xe(i))/d2(i));
    a12(i)= ((y1-ye(i))/d1(i))-((y2-ye(i))/d2(i));
    a13(i)= ((z1-ze(i))/d1(i))-((z2-ze(i))/d2(i));
    a21(i)= ((x1-xe(i))/d1(i))-((x3-xe(i))/d3(i));
    a22(i)= ((y1-ye(i))/d1(i))-((y3-ye(i))/d3(i));
    a23(i)= ((z1-ze(i))/d1(i))-((z3-ze(i))/d3(i));
    a31(i)= ((x1-xe(i))/d1(i))-((x4-xe(i))/d4(i));
    a32(i)= ((y1-ye(i))/d1(i))-((y4-ye(i))/d4(i));
    a33(i)= ((z1-ze(i))/d1(i))-((z4-ze(i))/d4(i));
    A= [a11(i) a12(i) a13(i);
        a21(i) a22(i) a23(i);
        a31(i) a32(i) a33(i)];
    f1v(i)= d2(i)-d1(i);
    f2v(i)= d3(i)-d1(i);
    f3v(i)= d4(i)-d1(i); 
    D= [(d21-f1v(i));
        (d31-f2v(i));
        (d41-f3v(i))];
    e= [eps21; eps31; eps41];
    R= [0.32932 -0.2018 -0.08995;
        -0.2018 0.3292 0.10616;
        -0.08995 0.10616 0.1143];
    DO= inv([(A.')*(inv(R))*A])*(A.')*(inv(R))*D;
    A*DO == D-e;
    dox(i)= DO(1,:);
    doy(i)= DO(2,:);
    doz(i)= DO(3,:);
    x(i)= xe(i)+dox(i);
    y(i)= xe(i)+doy(i);
    z(i)= ze(i)+doz(i);
    xe(i+1)= xe(i)+dox(i);
    ye(i+1)= xe(i)+doy(i);
    ze(i+1)= ze(i)+doz(i);
end
toc