>> global a b c d e v; a=1;b=0;c=-1;d=0;e=1.5;v=1.5; syms x y; [x,y]=solve('x^2+y^2-2*a*x-2*b*y=e^2-a^2-b^2','x^2+y^2-2*c*x-2*d*y=v^2-c^2-d^2'); x=vpa(x,4) y=vpa(y,4) x = (a^2 + b^2 - 1.0*c^2 - 1.0*d^2 - 1.0*e^2 + v^2 - (1.0*b*(a^2*b + b*c^2 + a^2*d - 1.0*b*d^2 - ... y = (0.5*(a^2*b + b*c^2 + a^2*d - 1.0*b*d^2 - 1.0*b^2*d - 1.0*b*e^2 + c^2*d + d*e^2 + b*v^2 - .... >> x=subs(x) x = 0 0 >> y=subs(y) y = -1.1180339887498948482045868343656 1.1180339887498948482045868343656 >> x=eval(x) x = 0 0 >> y=eval(y) y = -1.1180 1.1180 >>