El método estándar homogeneiza primero el problema, por ejemplo, estableciendo $y=x^3+u$ para que $u(0)=0=u(1)$ . La ecuación cambia a $$ 0=(6x+u'')+(3x+\tfrac1xu) +(x^3+u)^2=u''+\tfrac1xu+2x^3u+x^6+9x+u^2 $$
En los comentarios has identificado el operador lineal $L[u]=u''+\frac1xu'$ como el que hay que invertir por función de Green. El siguiente paso es determinar las soluciones de base de $L[u]=0\iff (xu')'=0\iff u=C\ln(x)+D$ que satisfacen una de las condiciones de contorno homogéneas, por ejemplo $u_1(0)=1$ , $u_1'(0)=0$ y $u_2(1)=0$ , $u_2'(1)=1$ . Esto da $u_1(x)=1$ y $u_2(x)=\ln(x)$ .
Entonces $\tilde G(s,x)=u_1(\min(s,x))u_2(\max(s,x))$ es continua y $\tilde G_x(s,x)=\chi_{[s,1]}(x)\frac1x$ tiene un salto de $0$ à $\frac1s$ en $x=s$ por lo que la función de Green se obtiene reduciendo ésta al salto unitario en $$G(s,x)=s\ln(\max(s,x)).$$ El operador asociado es $$G[u](x)=\int_0^1G(x,t)u(t)\,dt=x\ln(x)\int_0^xu(s)\,ds+\int_x^1s\ln(s)\,u(s)\,ds$$ para que $L[G[u]])=u$ . Aplicando esto a la ecuación diferencial completa se obtiene la ecuación integral de punto fijo $$ u=-G[2x^3u+x^6+9x+u^2] $$ Habría que comprobar si la singularidad en $x=0$ impide la finitud para las funciones continuas $u$ con $u(0)=u(1)=0$ .