Find the inverse of Interpolating function from NDSolve

Question

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.

Answer 1

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]

Answer 2

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 *)