A mundane but effective approach is to use the `Shooting` `Method` built into `NDSolve`.

    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`.

The `k = .4` solution, compared to that given in the question is, 

    LogPlot[{sol[x], solt[x]}, {x, x0, xMax}, AxesLabel -> {x, y}]

![enter image description here][1]

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 

    sol[x0]
    (* 4.51095 *)

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

![enter image description here][2]

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

![enter image description here][3]


  [1]: https://i.sstatic.net/oqqsL.png
  [2]: https://i.sstatic.net/Rs7j6.png
  [3]: https://i.sstatic.net/jmApx.png