Kernel quits while using ParametricNDSolve

Question

I'm trying to reproduce the plots (Figure 1 and Figure 2) of this paper. The plots are obtained by solving the Eqs. 4.3 and 4.4. The boundary conditions described in section 3.3 for Eqn. 4.3 are provided as follows: $t'(u_0)= \infty$, $t'(u_k)=\pm \infty$ where $u_0$ is the turning point of the surface and $u_k$ denotes the deep IR. Similarly, the boundary conditions for Eqn. 4.4 are provided as, $t'(\infty)= 0$, $t'(u_k)=\pm \infty$. For all the plots, they have fixed $u_k=1$.

I'm using this code for Eqn. 4.3,

     uk = 1;
eom1 = t'[u] == 
   Sqrt[(u0^2 uk^3 - 
     u0^5)/((u^3 - uk^3) (u0^2 uk^3 - u0^5 - u^2 uk^3 + u^5))];
sol = ParametricNDSolve[{eom1, t'[u0] == 10^3, 
   t'[uk] == 10^3}, {t}, {u, 0, 2}, {u0}]
Plot[Evaluate[Im[t[u0][1] /. sol]], {u0, 0, 2}]

but it's not working.

Can anybody suggest how to proceed with the code?

Answer 1

We can use NIntegrate instead of NDSolve as follows

uk = 1;(*b=u0^2 uk^3-u0^5*);
t[u1_?NumericQ, u2_?NumericQ, b_?NumericQ, t0_, pm_] := 
 pm  NIntegrate[
    Sqrt[(b)/((u^3 - uk^3)  (b - u^2  uk^3 + u^5))], {u, u2, uk, u1}, 
    Method -> "LocalAdaptive", Exclusions -> {u - uk == 0}] + t0

Then we can play with parameters and Re,Im to plot parts of trajectory, for example

plot1 = Quiet@
  ParametricPlot[{{u1, Re[t[u1, -2, 1, 0, 1]]}, {u1, 
     Re[t[u1, -2, 1, 0, -1]]}}, {u1, 1, 2}, PlotStyle -> {Blue, Blue},
    PlotRange -> {{0, 2}, {-2, 2}}, AspectRatio -> 1/2];

plot2 = Quiet@
  ParametricPlot[{{u1 - .19, Im[t[u1, 1, -1, -.5  I, 1]]}, {u1 - .19, 
     Im[t[u1, 1, -1, .5  I, -1]]}}, {u1, 1.19, 2.2}, 
   PlotStyle -> {Green, Green}, PlotRange -> {{0, 2}, {-2, 2}}];

Show[ParametricPlot[{uk, u}, {u, -2, 2}, PlotStyle -> {Black, Thin}, 
  PlotRange -> {{.8, 2}, {-2, 2}}, AspectRatio -> 1/2], plot1, plot2] 

Update 1. To compute red lines we use

uk = 1;
t[u1_?NumericQ, u2_?NumericQ, b_?NumericQ, t0_, pm_] := 
  pm NIntegrate[
     Sqrt[(b)/((u^3 - uk^3) (b - u^2 uk^3 + u^5))], {u, u2, uk, u1}, 
     Method -> "LocalAdaptive", Exclusions -> {u - uk == 0}] + t0;

plot = Table[
  Quiet@ParametricPlot[
    Table[{u1, 
      Re[t[u1, x, x^2 uk^3 - x^5, 0, m]]}, {m, {-1, 1}}], {u1, 1, x}, 
    PlotStyle -> Red, PlotRange -> {{0.5, 5}, {-2, 2}}, 
    AspectRatio -> 1/2], {x, 1.5, 4, .5}]

Also we can use NDSolve as follows

sol[u0_] := 
 Module[{uk = 1, e = 10^-5}, 
  s = NDSolve[{t'[
        u]^2 == (u0^2 uk^3 - 
         u0^5)/((u^3 - uk^3) (u0^2 uk^3 - u0^5 - u^2 uk^3 + u^5)), 
     t[u0 - e] == 0}, t, {u, uk + e, u0 - e}]; s]

Table[Quiet@
  ParametricPlot[
   Table[{u1, Re@Evaluate[t[u1] /. sol[x][[m]][[1]]]}, {m, 2}], {u1, 
    1, x}, PlotStyle -> Red, PlotRange -> {{0.5, 5}, {-2, 2}}, 
   AspectRatio -> 1/2], {x, 1.5, 4, .5}]