x0 = .01; xMax = 11.5; k = .4;
sol = NDSolveValue[{y''[x] + 3/x y'[x] - y[x] + 3/2 y[x]^2 - 
     k/2 y[x]^3 == 0, y'[xMax] == -y[xMax], y'[x0] == 0}, y, {x, x0, xMax}, 
  Method -> {"Shooting", "StartingInitialConditions" -> {y[x0] == 4.5, y'[x0] == 0}}]

Note that the boundary condition y'[xMax] == -y[xMax] is employed to match the asymptotic solution to the ODE, c Exp[-x]. Also, xMax is increased to 11.5.

Two features are evident. First, the solution already is exponentially decreasing at x = 3. Second, both approaches become unstable, the method in the question at about x = 9.5 and that given here at about x = 11. This presumably occurs due to unavoidable round-off errors that couple into the exponentially growing asymptotic solution to the ODE.

The value of y near the origin agrees quite closely with that in the question.

The next plot depicts y[x0] (points) as a function of k, along with the bounds on y[x0] (solid curves) given in the question. (These bounds correspond to the values of y[x0] which satisfy the ODE for all x, apart from the boundary condition at infinity.)

Solutions for k < .8 are obtained with moderate ease by first running the code above with xMax = 5, which is fairly forgiving of inaccurate initial guesses for y[x0], and then using the resulting actual y[x0] with xMax = 10 to obtain y[x0] accurate to five significant figures or better. The process becomes rather more delicate at larger k, because y[0] lies very close to its upper bound. I typically would need to guess about five times (moving the point at which NDSolve declared the equations stiff toward xMax)and also increase xMax somewhat to obtain a solution that was both stable and accurate. Any particular calculation takes less than a second of computer time.

x0 = .01; xMax = 11.5; k = .4;
c = N[D[(BesselK[1, x]/x), x]/(BesselK[1, x]/x) /. x -> xMax];
sol = NDSolveValue[{y''[x] + 3/x y'[x] - y[x] + 3/2 y[x]^2 - 
     k/2 y[x]^3 == 0, y'[xMax] == c y[xMax], y'[x0] == 0}, y, {x, x0, xMax}, 
  Method -> {"Shooting", "StartingInitialConditions" -> {y[x0] == 4.5, y'[x0] == 0}}]

Note that the boundary condition y'[xMax] == c y[xMax] is employed to match the exponentially decreasing asymptotic solution to the ODE, BesselK[1, x]/x. Also, xMax is increased to 11.5.

Two features are evident. First, the solution already is exponentially decreasing at x = 3. Second, both approaches become unstable, the method in the Question at about x = 9.5 and that given here at about x = 11. This presumably occurs due to unavoidable round-off errors that couple into the second, exponentially growing, asymptotic solution to the ODE.

The value of y near the origin agrees quite closely with that in the Question.

The next plot depicts y[x0] (points) as a function of k, along with the bounds on y[x0] (solid curves) given in the Question. (These bounds correspond to the values of y[x0] which satisfy the ODE for all x, apart from the boundary condition at infinity.)

data = {{0., 5.78018}, {0.1, 5.49104}, {0.2, 5.18459}, {0.3, 4.85871}, {0.4, 4.51095}, 
        {0.5, 4.13862}, {0.6, 3.73916}, {0.7, 3.31151}, {0.75, 3.08818}, 
        {0.8, 2.86076}, {0.85, 2.63322}, {0.88, 2.49923}, {0.9, 2.41187}};
Show[ListPlot[data, AxesLabel -> {"k", y[0]}, PlotStyle -> Black], 
  Plot[{(3 - Sqrt[9 - 8 k])/(2 k), (3 + Sqrt[9 - 8 k])/(2 k)}, {k, 0, 1},
  PlotStyle -> Black], PlotRange -> {{0, 1}, {0, 6}}]

Solutions for k < .8 are obtained with moderate ease by first running the code above with xMax = 5, which is fairly forgiving of inaccurate initial guesses for y[x0], and then using the resulting actual y[x0] with xMax = 10 to obtain y[x0] accurate to five significant figures or better. The process becomes rather more delicate at larger k, because y[0] lies very close to its upper bound. I typically would need to guess about five times (moving the point at which NDSolve declared the equations stiff toward xMax) and also increase xMax somewhat to obtain a solution that was both stable and accurate. Any particular calculation takes less than a second of computer time.

Two features are evident. First, the solution is exponentially decreasing beyond x = 3. Second, both approaches become unstable, the method in the question at about x = 9.5 and that given here at about x = 11. The value of y near the origin is

y[x0] can be computed (with xMax = 10) without much difficulty for k less than 0.8.

A comparison among y for k of 0(blue) and .6(brown), and .8(green) follows.

Two features are evident. First, the solution already is exponentially decreasing at x = 3. Second, both approaches become unstable, the method in the question at about x = 9.5 and that given here at about x = 11. This presumably occurs due to unavoidable round-off errors that couple into the exponentially growing asymptotic solution to the ODE.

The value of y near the origin agrees quite closely with that in the question.

The next plot depicts y[x0] (points) as a function of k, along with the bounds on y[x0] (solid curves) given in the question. (These bounds correspond to the values of y[x0] which satisfy the ODE for all x, apart from the boundary condition at infinity.)

Solutions for k < .8 are obtained with moderate ease by first running the code above with xMax = 5, which is fairly forgiving of inaccurate initial guesses for y[x0], and then using the resulting actual y[x0] with xMax = 10 to obtain y[x0] accurate to five significant figures or better. The process becomes rather more delicate at larger k, because y[0] lies very close to its upper bound. I typically would need to guess about five times (moving the point at which NDSolve declared the equations stiff toward xMax)and also increase xMax somewhat to obtain a solution that was both stable and accurate. Any particular calculation takes less than a second of computer time.

A comparison among y for k of 0(blue), .6(brown), .8(green), and .9 (red) follows. That the k = 0.9 solution is essentially constant out to about x = 10 explains why the calculation at larger k is so delicate and also why xMax = 15 is necessary for accuracy.

Note that I have not experimented with NDSolve Methods for handling stiffness, which may improve the robustness of this solution. Nonetheless, I was able to obtain the curves above without excessive difficulty.