Your function has a square root term of negative value:

In[182]:= Clear["Global`*"]

In[183]:= Simplify[
 y[x]^3 (1 - y[x]) x^2 - (1 - y[x])^3 x^2 + 3/8 y[x] (1/2 y[x]^(-1/2))]

Out[183]= 
x^2 (-1 + y[x])^3 + (3 Sqrt[y[x]])/16 - x^2 (-1 + y[x]) y[x]^3

In[184]:= s = 
 NDSolve[{y'[x] == 
    x^2 (-1 + y[x])^3 + (3 Sqrt[Abs[y[x]]])/16 - 
     x^2 (-1 + y[x]) y[x]^3, y[0] == 0}, y, {x, 0.000, 1}, 
  Method -> Automatic]

Out[184]= {{y -> InterpolatingFunction[{{0., 1.}}, <>]}}

It works if I take absolute value underneath the square root. I also used Automatic as the integration method.

May I make explicit a point which was implicit already in David's answer: The 'normal' way to obtain numerical solutions for ODEs from Mathematica is to call NDSolve and leave it to the system to find out among the many built-in methods and merthod variants the one which seems to be most suited for your equation. Since it definitively has methods which are designed to cope with stiff equations it should come up with a reliable solution whenever it exists mathematically. Of course, it can only make you learn and understand if you compare the behavior of standard methods for cases in which some of them get into trouble. Unfortunately the documentation on the 'Method' option of NDSolve is poor, so that some experimentation is needed to see which methods are implemented.