2D root finding algorithm using a recursive bisection method for convex quadrilaterals wich allows to find more than 1 zero

Lately I've been occupied in finding critical points in a 2D vector field. These are locations on a X, Y domain where both components of a vector field go to zero, i.e., (U, V) = (0,0). This is the same has having two functions, f(x,y) and g(x,y), and I want to find the set positions (xp,yp) where simultaneously f(xp,yp)=0 and g(xp,yp)=0.
This is not a topic I'm studying and I end up finding these points by mapping the domain with streamlines and verifying for patterns which indicate that a critical point is present (either saddles, nodes or focus). By verifying I meant visual inspection, or creating streamlines and looking at those, which is very boring.

Fortunately I've found the work of A. Globus, C. Levit, and T. Lasinski, where for the 2D case they do this process analytically by fitting bilinear interpolations and verifying the several particular cases (no solution inside the quadrilateral, 1 solution, 2 solutions, imaginary roots or infinite solutions).

Throughout this process I thought of a method of refining a root inside a quadrilateral, akin to the bisection method, but different. Instead of dividing to cut my quadrilateral in order to refine my solution, I create a tree where I further divide my element, letting me search for more roots in other locations. The best way to implement this was with a recursive function, as the quadrilateral divisions and root finding tests are the same to the children of the first element and subsequently. The test to see if a quadrilateral should be refined further is very basic: if the U values at the vertexes are all either >0 or <0, and the same happens to the V values at the vertexes, then this quadrilateral does not have any zero. Of course, there are functions where the topology would destroy this basic technique, such as one that has 4 roots inside the quadrilateral, located in a way that U and V is positive at all vertexes.

However for some functions this would allow to find more than one root. The vector field functions, U and V, may be discrete (leading to the bilinear interpolation case) or have an analytic representation. This would be achieved by changing the ufun and vfun.
Unfortunately it is not robust enough to sort out between repeated zeros (for example, if 2 or 4 children of a quadrilateral have the same vertex which coincides with a root, neither to work well under infinite solutions.

Either way it may be a useful piece of code for the future.