function [emf] = NewtonGPS4(deltaTi1)
format long;
% Define satellite positions - Part IV

% Position on Earth
x = 0;
y = 0;
z = 6370;
d = 0.0001;
u = [x; y; z; d];

% will change later
x0 = x;
y0 = y;
z0 = z;
d0 = d;

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

% Used for random positions; - code will change
i_1 = pi/3;
i_2 = pi/8;
i_3 = pi/5;
i_4 = pi/4;
i_5 = 3*pi/8;
i_6 = 3*pi/5;

phiMatrix = [i_1 i_2 i_3 i_4];
thetaMatrix = [i_1 i_2 i_5 i_6];
n = 4; %Use size of matrix;
conv = 5; % number of iteration to see convergence

rho = 26570;

phi_i = phiMatrix(1);
theta_i = thetaMatrix(1);

A1 = rho*cos(phi_i)*cos(theta_i);
B1 = rho*cos(phi_i)*sin(theta_i);
C1 = rho*sin(phi_i);   
R1 = sqrt(A1^2 + B1^2 + (C1 - 6370)^2);
t1 = d + R1/c;

phi_i =+ phi_i*0.05;
theta_i =+ theta_i*0.05;
A2 = rho*cos(phi_i)*cos(theta_i);
B2 = rho*cos(phi_i)*sin(theta_i);
C2 = rho*sin(phi_i);   
R2 = sqrt(A2^2 + B2^2 + (C2 - 6370)^2);
t2 = d + R2/c;

phi_i =+ phi_i*0.05;
theta_i =+ theta_i*0.05;
A3 = rho*cos(phi_i)*cos(theta_i);
B3 = rho*cos(phi_i)*sin(theta_i);
C3 = rho*sin(phi_i);   
R3 = sqrt(A3^2 + B3^2 + (C3 - 6370)^2);
t3 = d + R3/c;

phi_i =+ phi_i*0.05;
theta_i =+ theta_i*0.05;

A4 = rho*cos(phi_i)*cos(theta_i);
B4 = rho*cos(phi_i)*sin(theta_i);
C4 = rho*sin(phi_i);   
R4 = sqrt(A4^2 + B4^2 + (C4 - 6370)^2);
t4 = d + R4/c;

phi_i =+ phi_i*0.05;
theta_i =+ theta_i*0.05;
A5 = rho*cos(phi_i)*cos(theta_i);
B5 = rho*cos(phi_i)*sin(theta_i);
C5 = rho*sin(phi_i);   
R5 = sqrt(A5^2 + B5^2 + (C5 - 6370)^2);
t5 = d + R5/c;

phi_i =+ phi_i*0.05;
theta_i =+ theta_i*0.05;
A6 = rho*cos(phi_i)*cos(theta_i);
B6 = rho*cos(phi_i)*sin(theta_i);
C6 = rho*sin(phi_i);   
R6 = sqrt(A6^2 + B6^2 + (C6 - 6370)^2);
t6 = d + R6/c;

phi_i =+ phi_i*0.05;
theta_i =+ theta_i*0.05;
A7 = rho*cos(phi_i)*cos(theta_i);
B7 = rho*cos(phi_i)*sin(theta_i);
C7 = rho*sin(phi_i);   
R7 = sqrt(A7^2 + B7^2 + (C7 - 6370)^2);
t7 = d + R7/c;

phi_i =+ phi_i*0.05;
theta_i =+ theta_i*0.05;
A8 = rho*cos(phi_i)*cos(theta_i);
B8 = rho*cos(phi_i)*sin(theta_i);
C8 = rho*sin(phi_i);   
R8 = sqrt(A8^2 + B8^2 + (C8 - 6370)^2);
t8 = d + R8/c;

%error
%deltaTi1 = 10^-8;    
t1 = t1 - deltaTi1;
t2 = t2 + deltaTi1;
t3 = t3 - deltaTi1;
t4 = t4 + deltaTi1;
    
      for i = 1: conv   
       % 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
     
     u
     deltasPos = [x - u(1); y - u(2); z - u(3)];
     emf = norm(deltasPos, inf)/(c*(10^-8))
     maxEMF = max(emf)
end