c = 299792.458 ;
A1=15600;
B1=7540;
C1=20140;
A2=18760;
B2=2750;
C2=18610;
A3=17610;
B3=14630;
C3=13480;
A4=19170;
B4=610;
C4=18390;
t1=0.07074;
t2=0.07220;
t3=0.07690;
t4=0.07242;

ux = [(-2*A1+2*A2) ; (-2*A1+2*A3); (-2*A1+2*A4)];
uy = [(-2*B1+2*B2) ; (-2*B1+2*B3); (-2*B1+2*B4)];
uz = [(-2*C1+2*C2) ; (-2*C1+2*C3); (-2*C1+2*C4)];
ud = [-c^2*(-2*t1+2*t2); -c^2*(-2*t1+2*t3);-c^2*(-2*t1+2*t4)];
    
w = [((A1^2-A2^2)+(B1^2-B2^2)+(C1^2-C2^2)-c^2*(t1^2-t2^2));
((A1^2-A3^2)+(B1^2-B3^2)+(C1^2-C3^2)-c^2*(t1^2-t3^2));
((A1^2-A4^2)+(B1^2-B4^2)+(C1^2-C4^2)-c^2*(t1^2-t4^2))];
syms d;
x = (d*det([uy uz ud])+det([uy uz w]))/(-1*det([uy uz ux]));
y = (d*det([uz ux ud])+det([uz ux w]))/(-1*det([uz ux uy]));
z = (d*det([ux uy ud])+det([ux uy w]))/(-1*det([ux uy uz]));
f = (x-A1)^2+(y-B1)^2+(z-C1)^2-(c*(t1-d))^2;
format long
d = double(solve( f== 0, d))
x = (d*det([uy uz ud])+det([uy uz w]))/(-1*det([uy uz ux]));
y = (d*det([uz ux ud])+det([uz ux w]))/(-1*det([uz ux uy]));
z = (d*det([ux uy ud])+det([ux uy w]))/(-1*det([ux uy uz]));

% Here¡¯s an alternate form using Cramer¡¯s rule!
% ud = ud.*d
% D = det([ux uy uz]);
% Dx = det([(-w-ud) uy uz]);
% x= Dx/D;
% Dy = det([ux (-w-ud) uz]);
% y= Dy/D;
% Dz = det([ux uy (-w-ud)]);
% z= Dz/D;