Approximate inverse of a definite integral

I am still not sure what is the best (i.e. maximal accuracy) way to do that.
In this post @stevengj explains why it is better to approximate the function you want to integrate (if you are interested in a cumulative integration).

So I am trying that (following the example in the same thread).
I am also interested in the inverse function h^{-1}. So I use Newton’s Method to invert it for any value that I am interested in. I am not happy with the accuracy I get. Part of the complication, in my case, is that the function f(x) is actually a function that depends on another parameter f(x, a), and then for different a's I do not get consistent results.

I think the main problem comes from the approximation near the boundaries of the integral \int_{x_1}^{x_2} f(x') dx'. Is there any way to increase the approximation of ApproxFun.jl, or to account for the boundaries x_1, x_2 better?

Alternatively, is there a better way (more accurate) to compute this integral and its inverse?