Numerical bounce solution

The argument in the paper linked in the comments suggests that if one starts close enough to the one equilbrium at $y \approx c_1 = 1.04$, then the asymptotically vanishing solution may be found. Below is a solution that is not quite there. It starts at $x = 10^{-8}$ at $y = c_1 - 10^{-25}$ and uses WorkingPrecision -> 50; it certainly needs more than machine precision to ensure that $y$ does not round to $c_1$. It's not close enough. Getting close enough seems to require a very small difference and consequently a very high working precisions. But NDSolve fails when I push it, indicating an error test failure (NDSolveValue::nderr). It may be that for such an example to be computationally feasible, you need a greater difference in the value of V[x] at the minima.

ClearAll[V];
V[y_] = 1/4 (-1 + y)^2 y^2 - 25/10000 y^3;
ode = y''[x] + (3 y'[x])/x == V'[y[x]];

wp = 50;
cps = NSolve[V'[y] == 0, WorkingPrecision -> wp];

ndsol[x0_?NumericQ, x1_?NumericQ, y0_?NumericQ, yp_?NumericQ] := 
  NDSolveValue[{ode, y'[x1] == yp, y[x1] == y0, 
    WhenEvent[y[x] < 0, "StopIntegration"],
    WhenEvent[y[x] > 49/100 && x > x1, "StopIntegration"]}, 
   y, {x, x0, 10000}, 
   WorkingPrecision -> Precision@{x0, x1, y0, yp}];
ListLinePlot[
 ndsol[2^-50, 1/10^8, y - 10^-25 /. Last@cps, 0],
 PlotRange -> {All, {-1, 2}}]

---

Update

Here is a better result, FWIW. The method "Extrapolation" is more robust on this problem.

wp = 150;  (* working precision *)
cps = NSolve[V'[y] == 0, WorkingPrecision -> wp]; 
iObj[x0_?NumericQ, x1_?NumericQ, y0_?NumericQ, yp_?NumericQ] := 
 Catch@NDSolveValue[{ode, y'[x1] == yp, y[x1] == y0, 
    WhenEvent[y[x] < 0, Throw[x]],
    WhenEvent[y[x] > 49/100 && x > x1, Throw[$Failed]],
    WhenEvent[y[x] > 2 && x > x1, Throw[$Failed]]}, y, {x, x0, 1000}, 
   WorkingPrecision -> Precision@{x0, x1, y0, yp}, 
   Method -> "Extrapolation"];
obj[dy_?NumericQ] := iObj[2^-50, 1/10^8, y - 10^dy /. Last@cps, 0];

foo = Nest[  (* bisection method for finding optimal IC *)
  With[{mid = obj[Mean@#]},
    Switch[mid,
     $Failed, {First@#, Mean@#},
     _?NumericQ, {Mean@#, Last@#},
     _, {First@#, Mean@#}
     ]
    ] &,
  {-42 - 54/100, -42 - 53/100},
  60
  ]
(*
  {-(2452178951689226635963/57646075230342348800), 
   -(196174316135138130877/4611686018427387904)}
 *)
ndsol[x0_?NumericQ, x1_?NumericQ, y0_?NumericQ, yp_?NumericQ] := 
  NDSolveValue[{ode, y'[x1] == yp, y[x1] == y0, 
    WhenEvent[y[x] < 0, "StopIntegration"],
    WhenEvent[y[x] > 49/100 && x > x1, "StopIntegration"],
    WhenEvent[y[x] > 2 && x > x1, "StopIntegration"]}, 
   y, {x, x0, 10000}, WorkingPrecision -> Precision@{x0, x1, y0, yp}, 
   Method -> "Extrapolation"];
ListLinePlot[
 ndsol[2^-50, 1/10^8, y - 10^dy /. Last@cps, 0] /. dy -> First@foo,
 PlotRange -> All]

(Not an answer, an extended comment.)

So, I guess the method is not suitable for this problem. How can I get the numerical solution??"

Those statement and question do not make much sense to me. You already have the numerical solution. Why do you think the method is not suitable? The message InterpolationFunction::dmval is avoided if you decrease the range of tp a little (with 0.001). Study this result plot:

ListPlot[Table[{tp, Abs[Psol[tp, 50][50]]}, {tp, 0.8, 1.014, 0.001}], PlotRange -> All]

Also, maybe you should study the solution behavior for different values of the independent variable x not just 50.

For example, the following 3D plot gives an idea of the found solution behavior with respect to both tp and x. (We can see the "bouncing" along the x-axis.)

points = Flatten[Table[{tp, x, Abs[Psol[tp, 50][x]]}, {tp, 0.8, 1.014, 0.001}, {x, 1, 50, 1}], 1];
ListPlot3D[points, AxesLabel -> {"tp", "x", "y"}, PlotRange -> All]

OK, Michael E2's answer gives me the courage to post the following as an answer. It seems that finite difference method (FDM) can be used to find the solutions for this problem in a relatively easy way. I'll use pdetoae for the generation of finite difference equations:

Clear[V];
V[y_] = 1/4 (-1 + y)^2 y^2 - 25/10000 y^3;
ode = y''[x] + (3 y'[x])/x == V'[y[x]];    
newode = # x & /@ ode // Simplify
inf = 50;
domain = {0, inf};
bc = {y'[0] == 0, y[inf] == 0};
points = 500;
difforder = 2;
grid = Array[# &, points, domain];
(* Definition of pdetoae isn't included in this post,
   please find it in the link above. *)
ptoafunc = pdetoae[y[x], grid, difforder];
ae = ptoafunc[newode][[2 ;; -2]];
aebc = ptoafunc@bc;

initial[x_] := 2;
sollst = FindRoot[{ae, aebc}, Table[{y@x, initial@x}, {x, grid}]][[All, -1]];
ListLinePlot[sollst, DataRange -> domain]
sol = ListInterpolation[sollst, grid];
(* Error check: *)
Plot[Subtract @@ ode /. y -> sol // Evaluate, {x, 0, inf}, PlotRange -> All]

I hesitated to post this because:

1. Using different inf, I can find different solutions, e.g. inf = 100:

2. The error of the solution isn't that small, and even if I increase points, I cannot observe obvious error reduction:

   

But now, reading all the analysis above, I believe there just exist multiple solutions for the problem, and the code above is worth sharing.