# Finding solution for larger intervals with shooting method

I am trying to solve numerically the following non linear differential equation:

$$y''(x)+\frac{3}{x}y'(x)=\frac{{\rm d}V(y)}{{\rm d}y},\qquad V(y)=\frac{1}{4}(y(x)^2-1)^2+\frac{a}{2}(y(x)-1),$$

with boundary conditions:

$$y(\infty)= y_+,$$ $$\qquad$$ $$y'(0)=0,$$

where $$y_+$$ is the local minima of the potential $$V(y).$$

I am using the following code:

a = 0.23;
V[y_] := (1/4 )(y^2 - 1)^2 + a (y - 1)/(2);
ode = y''[x] + (3 y'[x])/x == V'[y[x]];
extrema = Block[{y}, y /. NSolve[D[V[y], y] == 0, y, Reals, WorkingPrecision -> 50]]

xs = 10^-8; xm = 18;
s = ParametricNDSolve[{ode, y'[xs] == 0, y[xm] ==extrema[[3]]}, y, {x, xs, xm}, {ys},
Method -> {"Shooting", "StartingInitialConditions" -> {y[xs] == ys, y'[xs] == 0}},
WorkingPrecision -> MachinePrecision,  AccuracyGoal -> 13, PrecisionGoal -> 13];

f = y[extrema[[1]] + 0.000001] /. s; (*choosing the optimal value of ys to get the correct solution*)
f[xs]
Plot[{f[x], extrema[[3]]}, {x, xs, xm}, AxesLabel -> {"x", "y(x)"},
PlotRange ->All]

With this code I manage to get a solution which satisfies my boundary conditions:

However, if I increase the interval, i.e. if I take xm larger (e.g. xm=20), then I am not able to get any solution and Mathematica gives the following error: