cyril - I think that there is a sign mismatch in the calculation of C. Rather than
C = y1-B/(y2-y1);
it should be
C = -y1+B/(y2-y1);
due to the re-arrangement of
x3(x2-x1) + y3(y2-y1) = B
=> y3 = [B - x3(x2-x1)]/(y2-y1)
=> y3 = B/(y2-y1) - x3(x2-x1)/(y2-y1)
and the above substitution into
(x3-x1)^2 + (y3-y1)^2 =r1^2
=> (x3-x1)^2 + (B/(y2-y1) - x3(x2-x1)/(y2-y1) - y1)^2 = r1^2
=> (x3-x1)^2 + (C - x3*tau)^2
where
C = B/(y2-y1) - y1
tau = (x2-x1)/(y2-y1)
Note that a couple of other things you could do is to replace the norm calls with just
r1 = (x3-x1)^2 + (y3-y1)^2; %norm(X(3,:) - X(1,:));
r2 = (x3-x2)^2 + (y3-y2)^2; %norm(X(3,:) - X(2,:));
to avoid the square root of norm and subsequent squaring (where needed) of r1 and r2. (Note that if you use the above two replacements, then you will have to remove the squares for r1 and r2.)
When I re-run your code with the above changes, I get the expected answer:
x =
-22.1513 -72.3327
x3 =
-22.1513
y =
34.5751 35.4920
y3 =
34.5751
