# Find the inverse of Interpolating function from NDSolve

I have a system of differential equations that I have numerically integrated to produce solutions. See below:

 ClearAll["Global*"]
dq = q'[z] == -(2 q[z]^2 + q[z] - 1)/(z + 1);
dy = y'[z] == (2 q[z]^2 + q[z] - 1)/(z + 1);
dh = h'[z] == h[z] (1 + q[z])/(z + 1);
dv =
v'[z] ==
(-x[z] (x[z] - q[z]) + 2 (x[z] + q[z]) - 3 (-v[z] + x[z] + q[z]) + 2 +
(-v[z] + x[z] + q[z]) (x[z] - 2 q[z] + 1) + 2 q[z]^2 + q[z] - 1)/(z + 1);
dx = x'[z] == -(-x[z] (x[z] - q[z]) + 2 (x[z] + q[z]) - 3 (-v[z] + x[z] + q[z]) + 2)/(z + 1);

sol = NDSolve[{dq, dy, dh, dv, dx, q == 0.499, y == 1 - 0.499, h == 54.0176,
v == 0.5, x == 0}, {q[z], y[z], h[z], v[z], x[z]}, {z, 20, 0}]


The output is the following: Plotting one of the solutions:

   R[z_] = (y[z] /. sol)
Plot[R[z], {z, 0, 20}]


Now I want to find the inverse of R[z] i.e. z[R]. I have tried InverseFunction:

  z[R_] = InverseFunction[R][z]


But this only gives me R^{-1}[z] as an output and when I try to plot z[R], nothing appears. Please help.

• The definition of R should read R[z_] = (y[z] /. sol[]); Then if you want to plot the inverse you can just use ParametricPlot[{R[z], z}, {z, 0, 20}, AspectRatio -> 1] Sep 12, 2022 at 22:44
• This question seems to be asked a lot: mathematica.stackexchange.com/… Sep 13, 2022 at 0:52

All of the follwing methods are almost equivalent.

• It is convience to use NDSolveValue and function names {q, y, h, v, x}.
Clear[sol, f, g];
sol = NDSolveValue[{dq, dy, dh, dv, dx, q == 0.499,
y == 1 - 0.499, h == 54.0176, v == 0.5,
x == 0}, {q, y, h, v, x}, {z, 20, 0}];
f = sol[];
g = InverseFunction[f];
Show[Plot[f[z], {z, 20, 0}, PlotRange -> All, PlotStyle -> Blue],
Plot[g[t], {t, f, f}, PlotRange -> All, PlotStyle -> Red],
AxesOrigin -> {0, 0}, AspectRatio -> Automatic]

• or use NDSolve and function names {q, y, h, v, x}.
Clear[sol, f, g];
sol = NDSolve[{dq, dy, dh, dv, dx, q == 0.499, y == 1 - 0.499,
h == 54.0176, v == 0.5, x == 0}, {q, y, h, v, x}, {z,
20, 0}];
f = y /. sol[]
g = InverseFunction[f];
Show[Plot[f[z], {z, 20, 0}, PlotRange -> All, PlotStyle -> Blue],
Plot[g[t], {t, f, f}, PlotRange -> All, PlotStyle -> Red],
AxesOrigin -> {0, 0}, AspectRatio -> Automatic]

• use the OP's setting. NDSolve and {q[z], y[z], h[z], v[z], x[z]} and we define f be the pure function by f = Function[z, y[z]] /. sol[].
Clear[sol, f, g];
sol = NDSolve[{dq, dy, dh, dv, dx, q == 0.499, y == 1 - 0.499,
h == 54.0176, v == 0.5, x == 0}, {q[z], y[z], h[z],
v[z], x[z]}, {z, 20, 0}]
f = Function[z, y[z]] /. sol[]
g = InverseFunction[f];
Show[Plot[f[z], {z, 20, 0}, PlotRange -> All, PlotStyle -> Blue],
Plot[g[t], {t, f, f}, PlotRange -> All, PlotStyle -> Red],
AxesOrigin -> {0, 0}, AspectRatio -> Automatic] The problem is that you need to extract the correct part of sol in some way.

For example,sol[[1, 1, 2, 0]] evaluates as InterpolatingFunction[...]. You can apply InverseFunction to this. For example

InverseFunction[sol[[1, 1, 2, 0]]][0.3]
(* 2.42446 *)
`