r = rootsb(f,g,xydomain)
finds the common zeros of two bivariate functions f and g in the domain xydomain (4-element vector), which are given as function handles.
This code exists besides roots(f,g) in chebfun2 (which does essentially the same task) because it is sometimes better for accuracy to resample the functions when working in a subdivided, smaller domain.
If xydomain is not provided, it defaults to [-1 1 -1 1].
test.m runs a simple test and shows the plots of the solution along with the zero curves. The solutions should be the intersections of the curves.
This code always employs the algorithm based on the hidden variable resultant method. For the algorithmic details, see
 Y. Nakatsukasa, V. Noferini, and A. Townsend, Computing the common zeros of two bivariate functions via Bezout resultants, submitted (2013).