# 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[20] == 0.499, y[20] == 1 - 0.499, h[20] == 54.0176,
v[20] == 0.5, x[20] == 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[[1]]); 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[20] == 0.499,
y[20] == 1 - 0.499, h[20] == 54.0176, v[20] == 0.5,
x[20] == 0}, {q, y, h, v, x}, {z, 20, 0}];
f = sol[[2]];
g = InverseFunction[f];
Show[Plot[f[z], {z, 20, 0}, PlotRange -> All, PlotStyle -> Blue],
Plot[g[t], {t, f[20], f[0]}, 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[20] == 0.499, y[20] == 1 - 0.499,
h[20] == 54.0176, v[20] == 0.5, x[20] == 0}, {q, y, h, v, x}, {z,
20, 0}];
f = y /. sol[[1]]
g = InverseFunction[f];
Show[Plot[f[z], {z, 20, 0}, PlotRange -> All, PlotStyle -> Blue],
Plot[g[t], {t, f[20], f[0]}, 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[[1]].
Clear[sol, f, g];
sol = NDSolve[{dq, dy, dh, dv, dx, q[20] == 0.499, y[20] == 1 - 0.499,
h[20] == 54.0176, v[20] == 0.5, x[20] == 0}, {q[z], y[z], h[z],
v[z], x[z]}, {z, 20, 0}]
f = Function[z, y[z]] /. sol[[1]]
g = InverseFunction[f];
Show[Plot[f[z], {z, 20, 0}, PlotRange -> All, PlotStyle -> Blue],
Plot[g[t], {t, f[20], f[0]}, 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 *)
`