I have an expression of the form:

expr = TerribleNumericalFunction[U[p]] == SimpleFunction[p] 

The Simple function is a polynomial, and the Terrible function involves elliptic integrals. I want to invert this, solving for U as a function of p as an interpolating function .

Naively, I would plug in p=0 and solve for U[0]. Then I could repeat this for many values of p and interpolate. Is there a better method? I feel like there should be a way to vectorize this process.

As a test case, this expression is way faster to evaluate than mine:

expr = Sqrt[0.04 + U[p] + U[p]^2] == 2.0 + 3.0*p; soln = Table[0, {i, 0, 99}]; For[p = 0, p < 100, p++, soln[[p]] = NSolve[expr, U[p]] ] 

This gets me a list of values (as substitution rules). This also gives 2 roots per p, which is true for my actual case.

How do I best turn my list of substitutions into a list of points for Interpolation? I guess there should really be two interpolating functions... one for each branch of U.