I have a system of two bivariate polynomials of degree 3 and would like to find its roots. With the help of resultants I project the system first at x-axis, then at y-axis, thus obtaining two univariate polynomial of degree 9 (Bezout's theorem). I know how to find the roots (or more precisely: the isolating intervals) of a univariate polynomial with arbitrary precision but I dont know how to make sure ("to certify") that a root of the first polynomial is also a root of the second polynomial.
How to find such pairs of intervals of these two univariate polynomials that will surely be also the roots of the original polynomial system?
First, you do not simply want to compare roots of univariate polynomials $f$ and $g$. For this task, you first determine isolating intervals for the roots of $f$ and $g$ and check whether they overlap. If no, the roots are distinct; otherwise compute $h = \operatorname{gcd}(f,g)$ and check the sign of (the square-free part of) $h$ at the endpoints of the intersection of the isolating intervals. If the signs are different, the corresponding roots of $f$ and $g$ are identical.
For certifying solutions of the bivariate system (say, $F = G = 0$), however, it will not suffice to compare roots of the resultant polynomials. In your example, you can determine the nine projections on the x-axis and nine on the y-axis; but this leaves you with a grid of 81 candidate solutions for the system. You will eventually be able to discard the non-solutions using interval arithmetic, if you refine the isolating intervals of the resultants sufficiently well: it will turn out that, for a non-solution, either $F \ne 0$ or $G \ne 0$. If you know a priori that the system is in generic position (that is, you have no solutions with higher multiplicity, such as tangential intersections between the algebraic curves described by the bivariate polynomials), and all solutions are real (or you work over the complex numbers), this could suffice.
Otherwise, you need to work harder to derive a certificate / stop condition / termination criterion for the actual solutions. This is non-trivial, and the key point of research on the Cylindrical Algebraic Decomposition (CAD) method. For the current state-of-the-art algorithm in that sector as well as references, have a look at On the complexity of computing with planar algebraic curves (disclaimer: I'm one of the authors).