A new numerical approach to compute all real roots of a system of two bivariate polynomial equations over a given box is described. Using the Bernstein-Bézier representation, we compute the best linear approximant and the best quadratic approximant of the two polynomials with respect to the L 2 norm. Using these approximations and bounds on the approximation errors, we obtain a fat line bounding the zero set first of the first polynomial and a fat conic bounding the zero set of the second polynomial. By intersecting these fat curves, which requires solely the solution of linear and quadratic equations, we derive a reduced subbox enclosing the roots. This algorithm is combined with splitting steps, in order to isolate the roots. It is applied iteratively until a certain accuracy is obtained. Using a suitable preprocessing step, we prove that the convergence rate is 3 for single roots. In addition, experimental results indicate that the convergence rate is superlinear (1.5) for double roots.