  We present a new algorithm for finding isolated zeros of a system of
real-valued functions in a bounded interval in $\mathbb{R}^n$. It uses the
Chebyshev proxy method combined with a mixture of subdivision, reduction
methods, and elimination checks that leverage special properties of Chebyshev
polynomials. We prove the method has R-quadratic convergence locally near
simple zeros of the system. We also analyze the temporal complexity and the
numerical stability of the algorithm and provide numerical evidence in
dimensions up to three that the method is both fast and accurate on a wide
range of problems. The algorithm should also work well in higher dimensions.
Our tests show that the algorithm outperforms other standard methods on this
problem of finding all real zeros in a bounded domain. Our Python
implementation of the algorithm is publicly available.
