I am looking to find all zeroes of a smooth function $ f \colon [0,1]^2 \to \mathbb{R}^2, $ using interval arithmetic.
The standard way to do this seems to be to use the interval Newton method
$$ x_{n+1} = x_n \cap \left ( y_n - \left ( \mathbf{J}f \right )^{-1}(x_n) f(y_n) \right ), $$
where $ y_n \in x_n $ (usually the midpoint).
This requires a way to compute and invert the interval Jacobian $\mathbf{J}f$. However, in my situation, the exact computation of the Jacobian is prohibitively expensive. Is there a variant on the interval Newton method which does not require the computation of the Jacobian? With normal arithmetic (instead of interval arithmetic), one has access to quasi-Newton methods such as Broyden's method, which estimate the inverse Jacobian as the method progresses, using function evaluations only.