Since you sought the zero contour, your implicit curve is presumably considered to be defined by the equation
(1-x)-a*x*y/(x^2+b*y^2) = 0.
You should not be afraid of solving for y in terms of x as a method of generating this curve, because doing so only involves solving a quadratic equation, whereas solving for x in terms of y would require solving a cubic equation with consequent greater complexity. The resulting quadratic equation will have two solutions
y1 = x*(a+sqrt(a^2-4*b*(1-x)^2))/(2*b*(1-x))
and (using the algebra trick of rationalizing the numerator)
y2 = 2*x*(1-x)/(a+sqrt(a^2-4*b*(1-x)^2))
both of which satisfy your equation with real numbers provided that
1-a/2/sqrt(b) <= x <= 1+a/2/sqrt(b).
As it turns out, the curve this equation defines is continuous throughout its entire infinite length, though it does cross over itself at the origin (0,0), and approaches the vertical line x = 1 asymptotically in both the y = inf and y = -inf directions.
It is not necessary to use either 'contour' or 'ezplot'. You can use the following code for this purpose. Note that to obtain these results it has been necessary to start with the y1 formula above with x moving forward, switch over to the y2 formula with x going backwards, and finally switch back to the y1 formula with x moving forward again, in order to trace a continuous path along the curve. This is forced by the nature of the curve and the fact that y is being solved in terms of x. If you look carefully at the curve, you will see why. (There is a white asterisk at the crossover point and the asymptote line x = 1 is shown in red.)
a = .1;
b = .001;
n = 3000;
ep = 1000*eps;
x1 = linspace(1.5,1+a/2/sqrt(b)-ep,n);
x2 = linspace(1+a/2/sqrt(b)-ep,1-a/2/sqrt(b)+ep,n);
x3 = linspace(1-a/2/sqrt(b)+ep,.5,n);
y1 = x1.*(a+sqrt(a^2-4*b*(1-x1).^2))./(2*b*(1-x1));
y2 = 2*x2.*(1-x2)./(a+sqrt(a^2-4*b*(1-x2).^2));
y3 = x3.*(a+sqrt(a^2-4*b*(1-x3).^2))./(2*b*(1-x3));
x = [x1,x2,x3];
y = [y1,y2,y3];
plot(x,y,'y-',0,0,'w*',[1,1],[y(1),y(end)],'r-')
You can check that the above satisfies the given equation throughout its length with the following:
f = (1-x)-a*x.*y./(x.^2+b*y.^2);
format long
max(abs(f))
