na.numerical analysis - 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.

• Next computer algebra - Basis for modular forms of half-integral weight
• Previous ag.algebraic geometry - Numerical evidence of Beilinson's conjecture in local fields and function fields