Finding the inverse of possibly non-invertible interpolation function

Question

Consider the following sample code

pts = Table[{i, i^2}, {i, -10, 10}];
foo = Interpolation[pts];
Plot[InverseFunction[foo][y], {y, 0, 10}]

From the plot, we see that InverseFunction[foo][y] only finds the positive x corresponding to the interpolation function foo, but ideally, I would like the inverse function to be something like ifoo[y]={x1,x2,...}, where x1,x2,... are the possible solutions for $y=x^2$.

In the end, I would like to generalize this to a multi-dimensional function. More specifically, let's say I have a set of data {{{x1,y1},{f1,g1}},...} and apply interpolation so that we get the function foo[x,y]={f[x,y],g[x,y]}. Now I would like to invert this function so that ifoo[z,w]={{x1,y1},{x2,y2},...} where {x1,y1},{x2,y2},... are the possible solutions towards foo[x,y]={z,w}.

How would one be able to do this?

Answer 1

Define the inverse with FindRoot

Clear["Global`*"]

pts = Table[{i, i^2}, {i, -10, 10}];
foo = Interpolation[pts];

inv[y_?NumericQ, branch_ : Automatic] :=
 Module[{x, init = If[branch === Automatic, -10, 10]},
  x /. FindRoot[foo[x] == y, {x, init}]]

Augmenting InverseFunction

Plot[{InverseFunction[foo][y], inv[y]},
 {y, 0, 100}, PlotLegends -> Placed["Expressions", {.85, .25}]]

Replacing InverseFunction

Plot[{inv[y, "+"], inv[y]}, {y, 0, 100}, 
 PlotLegends -> Placed["Expressions", {.85, .25}]]

Answer 2

You can find the two inverses by generating a differential equation of foo and solve with NDSolve applied to two initial conditions. Even derivatives of inverse function can be produced.

pts = Table[{i, i^2}, {i, -10, 10}];
foo = Interpolation[pts];

dinv = D[foo[x[y]] == y, y]

xsol[y_] = 
   x[y] /. NDSolve[{dinv, #}, x, {y, 100, 0}] & /@ {x[100] == -10, 
x[100] == 10} // Quiet

{Plot[xsol[y], {y, 0, 100}, ImageSize -> 300, AspectRatio -> 1], 
 Plot[xsol'[y], {y, 0, 100}, ImageSize -> 300, AspectRatio -> 1]
}