# NumericQ slows down the computation in ContourPlot

I am trying to plot the roots of a function f(x,y) using ContourPlot as follows

Sk[j_] = Sqrt[1 - x^2 + j 2 y + y^2];
ρ[j_] = (1 - x^2 + y^2 - Sk[j]^2)/(2 y);
A[s_, j_] = ( (s I x + Sk[j]) (y + ρ[j]))/(s - s y^2);
ζ[s_, j_] = ( (s I x + Sk[j]) (1 + y ρ[j]))/(s - s y^2);

Gs1[y_] = {{A[-1, 1], A[1, 1], A[-1, -1], A[1, -1]}, {1, 1, 1,
1}, {ζ[-1, 1], ζ[1, 1], ζ[-1, -1], ζ[
1, -1]}, {ρ[1], ρ[1], ρ[-1], ρ[-1]}};
Ms1[z_, y_] = {{E^(I Sk[1] z), 0, 0, 0}, {0, E^(-I Sk[1] z), 0,
0}, {0, 0, E^(I Sk[-1] z), 0}, {0, 0, 0, E^(-I Sk[-1] z)}};
Mxy = Inverse[Ms1[0, y]].Inverse[Gs1[y]].Gs1[y - 3].Ms1[0,
y - 3].Inverse[Ms1[8, y - 3]].Inverse[Gs1[y - 3]].Gs1[y].Ms1[8,
y];
ContourPlot[{N@Mxy[[4, 3]] == 0}, {x, -2, 2}, {y, 0.001, 3},
Exclusions -> y == 1, ExclusionsStyle -> None, PlotPoints -> 30,
ContourStyle -> {Directive[Blue]},
RegionFunction ->
Function[{x, y, f},  Abs[y - 1] - Abs[x] < -0.08 ]] // Timing


I had two problems:
1- In M10 (win10) it takes about 73.46 sec to get the results (see Fig. below) while in M12 (12.0.0.0 Linux x86 64-bit) I wait for more than 15 min but no results. So, why is this discrepancy?

2- To speed up the computation and enhance the results (smoothen the curves) I tried to define the function Mxy[[4, 3]] with numeric variables as follows

f[x_?NumericQ, y_?NumericQ] = Mxy[[4, 3]];

ContourPlot[{f[x, y] == 0}, {x, -2, 2}, {y, 0.001, 3},
Exclusions -> y == 1, ExclusionsStyle -> None, PlotPoints -> 30,
ContourStyle -> {Directive[Blue]},
RegionFunction ->
Function[{x, y, f},  Abs[y - 1] - Abs[x] < -0.08 ]] // Timing


but this takes a long time with the same output, see Fig. below! So, How can I speed up the computation and smooth the curves and why defining numeric variables did not help?

\$Version

(* "12.1.1 for Mac OS X x86 (64-bit) (June 19, 2020)" *)

Clear["Global*"]

Sk[j_] = Sqrt[1 - x^2 + j 2 y + y^2];
ρ[j_] = (1 - x^2 + y^2 - Sk[j]^2)/(2 y);
A[s_, j_] = ((s I x + Sk[j]) (y + ρ[j]))/(s - s y^2);
ζ[s_, j_] = ((s I x + Sk[j]) (1 + y ρ[j]))/(s - s y^2);

Gs1[y_] = {{A[-1, 1], A[1, 1], A[-1, -1], A[1, -1]}, {1, 1, 1,
1}, {ζ[-1, 1], ζ[1, 1], ζ[-1, -1], ζ[
1, -1]}, {ρ[1], ρ[1], ρ[-1], ρ[-1]}};
Ms1[z_, y_] = {{E^(I Sk[1] z), 0, 0, 0}, {0, E^(-I Sk[1] z), 0, 0}, {0, 0,
E^(I Sk[-1] z), 0}, {0, 0, 0, E^(-I Sk[-1] z)}};
Mxy = Inverse[Ms1[0, y]].Inverse[Gs1[y]].Gs1[y - 3].Ms1[0, y - 3].Inverse[
Ms1[8, y - 3]].Inverse[Gs1[y - 3]].Gs1[y].Ms1[8, y];


Simplify the expression prior to using it in ContourPlot

expr = Mxy[[4, 3]] // Simplify

(* -((3 E^(-8 I (Sqrt[-x^2 + (-4 + y)^2] - Sqrt[-x^2 + (-1 + y)^2])) (-1 + E^(
16 I Sqrt[-x^2 + (-4 + y)^2])) x (x +
I Sqrt[-x^2 + (-1 + y)^2]))/(2 Sqrt[-x^2 + (-4 + y)^2]
Sqrt[-x^2 + (-1 + y)^2] (-1 + y))) *)

ContourPlot[expr == 0, {x, -2, 2}, {y, 0.001, 3}, Exclusions -> y == 1,
ExclusionsStyle -> None, PlotPoints -> 30,
ContourStyle -> {Directive[Blue]},
RegionFunction ->
Function[{x, y, f}, Abs[y - 1] - Abs[x] < -0.08]] // Timing
`