Most efficient root finding algorithm for a monotonic function This is my first time asking a question here, so I may not be asking this in the right place. I am trying to find the roots of a monotonic function with as few function evaluations as possible. 
I have approximated a manifold with a piece-wise defined polynomial. The manifold is periodic and so I am only considering its unit cell (one period). I split the domain of the manifold, a parallelogram, into triangles. I then approximate each sheet of the manifold within each triangle with a unique quadratic polynomial. Here is the approximated manifold. I would like to find the root that satisfies the equation
\begin{equation}
\sum_{i}^{\mathrm{sheets}} \sum_{j}^{\mathrm{triangles}} \int p_{i,j}(x,y) \, \mathrm{d}C_{i,j} - A = 0
\end{equation}
where $i$ is a sum over the sheets of the manifold, $j$ is a sum over the triangular tiles, $p_{i,j}$ is the second degree polynomial approximation of the manifold's $i$th sheet within the $j$th tile, and $C_{i,j}$ is the region within a level curve of the polynomial approximation of the manifold's $i$th sheet within the $j$th tile. Here is a plot of the $C_{i,j}$ for each triangle and sheet for some estimate of the root. Said another way, I would like to find an isovalue where the area within the level curves of the polynomials, regions where the polynomials are less than the isovalue, is some predetermined value $A$.
At the moment I am using the bisection method, which is very slow because at each iteration it takes a significant amount of time to interpolate the manifold and then calculate the level curves and their containing areas. I may have hundreds of triangles and tens of sheets. I also tried the regular falsi method but ran into cases where its convergence was worse than the bisection method.
I was thinking that a bracketing method would work best. I also thought that I could take advantage of the fact that the function is monotonic.
 A: Regula falsi is known to have many drawbacks. As suggested by Lutz Lehmann, many improvements to regula falsi are known which recover at least some of its asymptotic order of convergence. One of the simplest is the Illinois method, but this does not recover the speed of the secant method and also fails to account for potentially bad behavior on a more global scale. Hybrid methods using 3-point interpolation are capable of handling these issues, but are more complicated to code. Some such methods include Brent's method, Dekker's method, and Chandrupatla's method. Of these, I recommend trying Chandrupatla's method.
Overall, not much can be taken advantage of from the fact that the searched function is monotone other than the existence of a unique root which can be bracketed. In particular, a crude initial interval can be found by simply searching further out (in the direction of $-f(x)$) until a sign change is detected. With a good hybrid method, extremely poor initial intervals are not a major issue, since they will close quickly with bisection until better methods may be used instead.
