The analysis of equilibria of ordinary differential equations (ODEs) that represent biochemical reaction networks is crucial in order to understand various functional properties of regulation in systems biology. In this paper, we develop a numerical algorithm to compute equilibria under the assumption that the regulatory network satisfies certain graph-theoretic conditions which lead to fixed-point iterations over an anti-monotonic function. Unlike generic approaches based on Newton’s method, our algorithm does not require the availability of the Jacobian of the ODE vector field, which may be expensive when the dimensionality of the system is large. More important, it produces an estimation (through over-approximation) of the entire set of equilibria, with the guarantee of yielding the unique equilibrium of the ODE in the case that the returned set is a singleton. We demonstrate the applicability of our algorithm to two signaling pathways of MAPK and EGFR.