In[]:=
U[y_,b_]:=-
2
b
-
2
b
2
y

(1+4
6
b
2
y
-6
4
y
+4
2
b
2
y
(-2+
4
y
)+
4
b
(-2+8
4
y
))
3/2
(
2
b
+
2
y
)
In[]:=
s=NDSolve[{f'[y](-U[y,1.6])^(1/2),f[-0.47]0},f,{y,-3,-0.47}]
Out[]=
fInterpolatingFunction
Domain: {{-3.,-0.47}}
Output: scalar

In[]:=
Plot[Evaluate[f[y]/.s],{y,-3,-0.47},PlotRangeAll]
Out[]=
In[]:=
V[y_,b_]:=-
4
-2
2
b
2
y

2
y
2
(-1+2
4
b
+2
2
b
2
y
)
2
b
+
2
y
+
2
b
-
2
b
2
y

(1+4
6
b
2
y
-6
4
y
+4
2
b
2
y
(-2+
4
y
)+
4
b
(-2+8
4
y
))
3/2
(
2
b
+
2
y
)
In[]:=
Plot[V[y,1.6],{y,-3,-0.47},PlotRange{{-3,-0.47},{-19,1.5}}]
Out[]=
In[]:=
ParametricPlot[{Evaluate[f[y]/.s],V[y,1.6]},{y,-3,-0.47},PlotRange{{-1.5,0},{-19,1}}]
Out[]=
-1.4
-15
-10
-5
In[]:=
ParametricPlot[{Evaluate[f[y]/.s],V[y,1.6]},{y,-3,-0.47},PlotRangeAll]
Out[]=
-1.0
-0.5
0.5
1.0
-1.0
-0.5
0.5
1.0
​