Updating Wagon's FindAllCrossings2D[] function

Here is my latest code for this function, from Chapter 12 of the third edition of "Mathematica in Action". It is pretty short, but I will let you work out if it is faster or more robust than yours. Note the PlotPoints option for difficult cases.

FindRoots2D::usage = 
  "FindRoots2D[funcs,{x,a,b},{y,c,d}] finds all nontangential solutions to
   {f=0, g=0} in the given rectangle."; 

Options[FindRoots2D] = {PlotPoints -> Automatic, MaxRecursion -> Automatic}; 

FindRoots2D[funcs_, {x_, a_, b_}, {y_, c_, d_}, opts___] := Module[
  {fZero, seeds, signs, fy}, 
  fy = Compile[{x, y}, Evaluate[funcs[[2]]]]; 

  fZero = Cases[Normal[
     ContourPlot[
        funcs[[1]] == 0, 
        {x, a-(b-a)/97, b+(b-a)/103}, {y, c-(d-c)/98, d+(d-c)/102}, 
        Evaluate[FilterRules[{opts}, Options[ContourPlot]]]]], 
     Line[z_] :> z, Infinity]; 

  seeds = Flatten[(
     (signs = Sign[Apply[fy, #1, {1}]]; 
      #1[[1 + Flatten[Position[Rest[signs*RotateRight[signs]], -1]]]]) &
     ) /@ fZero, 1];
  If[seeds == {}, {}, 
     Select[
        Union[({x, y} /.
           FindRoot[{funcs[[1]], funcs[[2]]}, {x, #1[[1]]}, {y, #1[[2]]}, 
              Evaluate[FilterRules[{opts}, Options[FindRoot]]]] & ) /@ seeds, 
           SameTest -> (Norm[#1 - #2] < 10^(-6) & )], 
        a <= #1[[1]] <= b && c <= #1[[2]] <= d & ]]]

This is ContourPlot based but seems much shorter:

FindCrossings2D[{f_, g_}, {x_, xmin_, xmax_}, {y_, ymin_, ymax_}] := 
  {x, y} /. (FindRoot[{f[x, y] == 0, g[x, y] == 0}, {{x, #[[1]]}, 
  {y, #[[2]]}}] & /@ (ContourPlot[{f[x, y] == 0, g[x, y] == 0}, 
  {x, xmin, xmax}, {y, ymin, ymax}][[1, 1]]))

It works:

f[x_, y_] := -Cos[y] + 2 y Cos[y^2] Cos[2 x];
g[x_, y_] := -Sin[x] + 2 Sin[y^2] Sin[2 x];

pts = FindCrossings2D[{f, g}, {x, -7/2, 4}, {y, -9/5, 21/5}];

ContourPlot[{f[x, y] == 0, g[x, y] == 0}, {x, -7/2, 4}, {y, -9/5, 
  21/5}, Epilog -> {AbsolutePointSize[6], Red, Point /@ pts}]

Let me give a different approach. FindRoot does a good job, but maybe we can calculate the seed-points in a different way. When you want to find the common roots of $f(x,y)$ and $g(x,y)$ you can transform the problem into one equation which has the same roots $$0=f(x,y)^2+g(x,y)^2$$ The nice property here, which I will use is that the right hand side of this equation is always greater than zero. The bad thing is, that functions that don't cross zero are kind of difficult for some numerical methods, but that will be no concern here.

When we look at our function and think of it as a kind of water-pool with a very low level of water, you would get something like this

Now the idea is to go around each of those small pools which have maybe one, maybe more local zeroes in it, and start from every coast point a root-search.

That's the time where the image-processing kicks in. Our function is always positive which gives a really nice image (I inverted gray-levels):

Cutting off the image at sea-level is just a binarization of the image.

Finding the coast-line of each pool is simply implemented by an image subtraction and a dilation of the binarized image.

The complete method is therefore to raster the above function, extract all coast-pixel with image processing and run FindRoot for each coast-point.

FindCrossings2D[{f_, g_}, {x_, xmin_, xmax_}, {y_, ymin_, ymax_}, 
                 n_, threshold_] := Module[
{seeds = ImageData[ImageSubtract[Dilation[#, 1], #] &@
  Binarize[ColorNegate[Image[
         Table[f[x, y]^2 + g[x, y]^2, 
          {y, ymin, ymax, (ymax - ymin)/(n - 1.0)}, 
          {x, xmin, xmax, (xmax - xmin)/(n - 1.0)}]]], 
       threshold], "Bit"]},
  DeleteDuplicates[Last@Last@Reap[MapIndexed[
    If[#1 === 1, Sow[{x, y} /. FindRoot[{f[x, y] == 0, g[x, y] == 0}, 
     {{x, Rescale[#2[[2]], {1, n}, {xmin, xmax}]},
      {y, Rescale[#2[[1]], {1, n}, {ymin, ymax}]}}]]] &, seeds, 2]
   ], (Norm[#1 - #2] < 10.^(-6)) &]
] 

Here n is the raster-size and thresh is the binarization threshold which should be a bit smaller than 1.

f[x_, y_] := -Cos[y] + 2 y Cos[y^2] Cos[2 x];
g[x_, y_] := -Sin[x] + 2 Sin[y^2] Sin[2 x];

roots = FindCrossings2D[{f, g}, {x, -7/2, 4}, {y, -9/5, 21/5}, 400, 0.8];

This approach has clearly the disadvantage of having a fixed raster size, while ContourPlot uses adaptive sampling. Nevertheless, for raster-sizes from 200-500 and thresholds from ?-0.95 the method finds all or at least many roots.