function [ u ] = NewtonGPS( )
% Position on Earth
x = 0;
y = 0;
z = 6370;
d = 0;
u = [x; y; z; d];

% First Trial #1
n = 100;

% Simultaneous (Ai, Bi, Ci) positions of satellites in km at ti time
A1 = 15600;
B1 = 7540;
C1 = 20140;
t1 = 0.07074;

A2 = 18760;
B2 = 2750;
C2 = 18610;
t2 = 0.07220;

A3 = 17610;
B3 = 14630;
C3 = 13480;
t3 = 0.07690;

A4 = 19170;
B4 = 610;
C4 = 18390;
t4 = 0.07242;

%Speed of light measured in km/sec
c = 299792.458;

x0 = x;
y0 = y;
z0 = z;
d0 = d;

    for i = 1: n
       % Find F matrix (4x1)
       F = [(x0 - A1)^2 + (y0 - B1)^2 + (z0 - C1)^2 - (c*(t1 - d0))^2;
           (x0 - A2)^2 + (y0 - B2)^2 + (z0 - C2)^2 - (c*(t2 - d0))^2;
           (x0 - A3)^2 + (y0 - B3)^2 + (z0 - C3)^2 - (c*(t3 - d0))^2;
           (x0 - A4)^2 + (y0 - B4)^2 + (z0 - C4)^2 - (c*(t4 - d0))^2];     
       % Find dF matrix (4x4)
       dF = [2*(x0 - A1)  2*(y0 - B1)  2*(z0 - C1)  2*(c^2*(t1 - d0));
       2*(x0 - A2)  2*(y0 - B2)  2*(z0 - C2) 2*(c^2*(t2 - d0));
       2*(x0 - A3)  2*(y0 - B3)  2*(z0 - C3)  2*(c^2*(t3 - d0));
       2*(x0 - A4)  2*(y0 - B4)  2*(z0 - C4)  2*(c^2*(t4 - d0))]; 
       % Newton Method to find next position
       v = dF\-F;
       u = u + v; 
       x0 = u(1);
       y0 = u(2);
       z0 = u(3);
       d0 = u(4);
    end
end

Supplier