Graphing Intersection of Two Implicit Functions when Equal to Zero

I have created a code that will use ezplot to graph the real and imaginary part of a wavefunction. I am interested in highlighting the intersection of the two parts (real and imaginary) when equal to zero. I have something similar, but I cannot figure out how to make it work only when the intersection is zero. I am using a scatter plot to overlay the red circles (the intersections which represent vortices) onto the graph. Currently, I have too many red circles. I have included the code below and have done a significant amount of online research for this, but I am wondering if there is not an easy way to do this.

I also have one small issue - I know ezplot has trouble graphing when there is no independent variable. In the case r = 1 below, the real and imaginary portions are:

psi_real0 = x^2+y^2-x psi_imag0 = -y

The real portion is being graphed correctly, but the imaginary part should be a horizontal line. I have tried adding a dummy variable such that psi_imag0 = 0*x-y, but this doesn't seem to be working.

The intersection is more important than this last part, but any help would be great.

By the way, this forum has been very helpful for this project, thanks in advance for all of the input.

Thank You, Ima

clc;
clear;
syms x y t k j e r;

outtype = 1; %1 for snapshots, 2 for movies
timescale = 15

N = 10
e = 1
psi_0old = expand((x + j*y)*(x - e - j*y))
psi_0 = subs(psi_0old, [j^2,j^3,j^4,j^5,j^6,j^7,j^8], [-1,-j,1,j,-1,-j,1])

lap_psi_0 = laplacian(psi_0)

psi_old = expand((t*j/2*lap_psi_0)+psi_0)

psi_real0 = 0*x+x^2+y^2-x+e^6*(t*(- x^4*y - x^2*y + 2*x*y^3 - y^5 - y^3 - 2*x^2*y^3 + 2*x^3*y)+t^2*(-8*x^4 + 16*x^3 - 16*x^2*y^2 - 10*x^2 + 16*x*y^2 + 2*x - 8*y^4 - 2*y^2)+t^3*(20*y/3+t^4*(12)))
psi_imag0 = 0*x+(2*t-y)+e^6*(t*(x^6 - 3*x^5 + 3*x^4*y^2 + 3*x^4 - 6*x^3*y^2 - x^3 + 3*x^2*y^4 + 4*x^2*y^2 - 3*x*y^4 - x*y^2 + y^6 + y^4)+t^2*(8*x*y - 8*x^2*y - 8*y^3 - 2*y)+t^3*(-20*x^2 + 20*x - 20*y^2 - 4))

if outtype==2
    for r = 1:N+1 
    psi_real=subs(collect(psi_real0,t),t,(r-1)/timescale)
    p = ezplot(psi_real)
    set(p, 'Color', 'g');
    hold on;
    psi_imag=subs(collect(psi_imag0,t),t,(r-1)/timescale)
    h = ezplot(psi_imag)
    set(h, 'Color', 'b');
    psi_diff = psi_real-psi_imag;
    [x,y]=solve(psi_real,psi_imag,x,y)
    scatter(x,y,'r')
    hold off;
    grid on;
    title(['t = ' num2str((r-1)/timescale)]);
    Z(r) = getframe(gcf);
    end
else

for r = 1:N+1 
   switch r
     case 1   
     subplot(2,2,1)    
    dr=1;
    case 2
        subplot(2,2,2)    
    dr=1;
      case 3
        subplot(2,2,3)    
    dr=1;
      case 4
        subplot(2,2,4)    
    dr=1;
       otherwise
    dr=0;
     end

if dr==1

psi_real=subs(collect(psi_real0,t),t,(r-1)/timescale)
p = ezplot(psi_real)
set(p, 'Color', 'g');
hold on;
psi_imag=subs(collect(psi_imag0,t),t,(r-1)/timescale)
h = ezplot(psi_imag)
set(h, 'Color', 'b');
[x,y]=solve(psi_real,psi_imag)
scatter(x,y,'r')
hold off;
grid on;
axis([-3 3 -3 3])
title(['t = ' num2str((r-1)/timescale)]);
Z(r) = getframe(gcf);
dr==0;
end

end

end

[ax4,h3]=suplabel(['Two Vortices: Non-Linear (+)(-)']  ,'t');
set(h3,'FontSize',15)