syms x y
a = 3;
b = 4;
c = 1;
d = 2;
g = 2.5;
eqn = (sqrt((x-a).^2 + (y-b).^2) - sqrt((x-c).^2 + (y-d).^2)).^2 == g.^2;
sol = solve(eqn,y);
xnum = -20:0.1:20;
ynum = double(subs(sol,x,xnum));
figure(1)
plot(xnum,real([ynum(1,:);ynum(2,:)]))
figure(2)
plot(xnum,imag([ynum(1,:);ynum(2,:)]))