function position = SolveGPS(satellites)
A = zeros(3, 3);
b = zeros(3, 2);
speed_of_light = 299792.458;
speed_of_light_squared = speed_of_light^2;

for i = 1 : 3
    left_constants = 0;
    for j = 1 : 3
        A(i, j) = satellites(1, j) - satellites(i + 1, j);
        left_constants = left_constants + (satellites(i + 1, j)^2 - satellites(1, j)^2);
    end
    left_constants = 0.5*left_constants;
    b(i, 1) = speed_of_light_squared*(satellites(1, 4) - satellites(i + 1, 4));
    b(i, 2) = 0.5*speed_of_light_squared*(satellites(i + 1, 4)^2 - satellites(1, 4)^2) - left_constants;
end

x = A\b;

positions = SolveGPSQuadratic(satellites, x);
if norm(positions(1, 1 : 3)) < norm(positions(2, 1 : 3))
    position = positions(1, :);
else
    position = positions(2, :);
end

