[x1,y1,x2,y2]=solve(\'x1^2 + y1^2=1\',\'x2^2-8*x2 +y2^2 +15=0\',\'x1*x2 + y1 * y2=1\',\'x1 + x2 =a\',\'x1\',\'y1\',\'x2\',\'y2\')

[x1,y1,x2,y2] = solve(...
x1^2/r1^2 + y1^2/r2^2 == 1,...
(x2-a)^2 + (y2-b)^2 == r3^2,...
x1*x2/r1^2 + y1*y2/r2^2 == 1,...
(x2-x1)^2 + (y2-y1)^2 + (x2-a)^2 + (y2-b)^2 == (x1-a)^2 + (y1-b)^2,...
x1,y1,x2,y2)

cond1 = a > 0
cond2 = b > 0
cond3 = r1 > 0
cond4 = r2 > 0
cond5 = r3 > 0
 
syms x1 y1 x2 y2 r1 r2 r3 a b 
eqn = [ x1^2 + y1^2 == r1^2 , (a - x2)^2 + (b - y2)^2 == r3^2, (x1*x2) + (y1*y2) == r1^2, (x1 - x2)^2 + (y1 - y2)^2 + (a - x2)^2 + (b - y2)^2 == (a - x1)^2 + (b - y1)^2]
sol = solve(eqn, [x1, y1, x2, y2])

eqn =
 
[ x1^2/r1^2 + y1^2/r2^2 == 1, (a - x2)^2 + (b - y2)^2 == r3^2, (x1*x2)/r1^2 + (y1*y2)/r2^2 == 1, (x1 - x2)^2 + (y1 - y2)^2 + (a - x2)^2 + (b - y2)^2 == (a - x1)^2 + (b - y1)^2]