Multivariate Bisection cross post in StackOverflow
I need an algorithm to perform a 2D bisection method for solving a $2$x$2$ non-linear problem. Example: two equations $f(x,y)=0$ and $g(x,y)=0$ which I want to solve simultaneously. I have very familiar with the 1D bisection ( as well as other numerical methods ). Assume I already know the solution lies between the bounds $x_1 < x < x_2$ and $y_1 < y < y_2$.
In a grid the starting bounds are:
    ^
    |   C       D
y2 -+  o-------o
    |  |       |
    |  |       |
    |  |       |
y1 -+  o-------o
    |   A       B
    o--+------+---->
       x1     x2

and I know the values at $f(A)$, $f(B)$, $f(C)$ and $f(D)$ as well as $g(A)$, $g(B)$, $g(C)$ and $g(D)$. I might even know for which edges $f=0$ and for which $g=0$.
To start the bisection I guess we need to divide the points out along the edges as well as the middle.
    ^
    |   C   F   D
y2 -+  o---o---o
    |  |       |
    |G o   o M o H
    |  |       |
y1 -+  o---o---o
    |   A   E   B
    o--+------+---->
       x1     x2

Now considering the possibilities of combinations such as checking if $f(G)*f(M)<0$ AND $g(G)*g(M)<0$ seems overwhelming. Maybe I am making this a little too complicated, but I think there should be a multidimensional version of the Bisection, just as Newton-Raphson can be easily be multidimed using gradient operators.
Any clues, comments, or links are welcomed.
 A: A remarkable generalization of bisection to multiple dimensions is the subgradient method from convex optimization theory. If $f$ and $g$ are convex then $h = f^2 + g^2$ is also convex, and a simultaneous zero of $f$ and $g$ is a minimum of $h$.
Unfortunately, the subgradient method has more theoretical than practical value. But in a two-dimensional problem, it might do okay.
A: You might want to consider the vector field
$ \vec{F}(x,y) = (f(x,y), g(x,y)) $
and look for sources and sinks of $\vec{F}$. I think this could be done by recursively dividing up the plane into squares and calculating the winding number of each square. If it is nonzero then you have a critical point within that square (cf Thm 2 this paper) and should divide further.
A: *

*Check the pair of opposite corners to determine if zeroes lie within each of the four subdivided rectangles (zeroes can be there in more than one of them). Eg. if f(M)>0 and f(A)<0, then AEMG contains zeroes of f. Same is true also if f(G)>0 and f(E)<0.

*Do this for all the four sub rectangles, and for both f and g.

*There will be atleast one which contains zeroes for both f and g. Zoom into that and repeat.
