Mathematica draw implicit function common tangent

implicit function common tangent.wls

(*
written in 2018/8/29
*)

Clear["Global`*"];
(*f[x,y]=0,(x1,y1),g[x,y]=0,(x2,y2) 本代码略去了两条曲线正好相切于同一点处的切线*)

f[x_, y_] = (x)^2 + x y + y^2 - 1;
g[x_, y_] = (x - 4)^2 + x y + (y)^2 - 1;(*两个隐函数*)

k1[x1_, y1_] = -(D[f[x1, y1], x1]/D[f[x1, y1], y1]);
k2[x2_, y2_] = -(D[g[x2, y2], x2]/D[g[x2, y2], y2]);
m1[x1_, y1_] = -(D[f[x1, y1], y1]/D[f[x1, y1], x1]);
m2[x2_, y2_] = -(D[g[x2, y2], y2]/D[g[x2, y2], x2]);
(*各自在 (x1,y1),(x2,y2) 处的切线斜率*)

a = NSolve[{f[x1, y1] == 0, g[x2, y2] == 0, 
    k1[x1, y1] == k2[x2, y2], 
        y2 == k1[x1, y1] (x2 - x1) + y1, 
        k2[x2, y2] (x1 - x2) + y2 == y1}, {x1, y1, x2, y2}, Reals];
b = NSolve[{f[x1, y1] == 0, g[x2, y2] == 0, 
    m1[x1, y1] == m2[x2, y2], 
        x2 == m1[x1, y1] (y2 - y1) + x1, 
        m2[x2, y2] (y1 - y2) + x2 == x1}, {x1, y1, x2, y2}, Reals];
(*先保证点在曲线上，然后确保切线是同一条*)
c = Union[a, b];
h[x_, y_] = ((y1*x2*(-1)) + (y2*x1) + (((x1*(-1)) + x2)*
            y) + ((y1 + (y2*(-1)))*x)) /. c;(*把解得的值代入*)c // Length
Show[ContourPlot[
    Evaluate@{f[x, y] == 0, g[x, y] == 0, h[x, y] == 0}, {x, -5, 
      5}, {y, -5, 5}, PlotLegends -> "Expressions"], 
  ListPlot[{{x1, y1} /. c, {x2, y2} /. c}]](*绘制切点*)