How to choose $y(0)$ if $y'(x)=-\alpha y^2(x)/x^2$? Say we have a differential equation $\frac{dy}{dx}=-\alpha \frac{y^2}{x^2}$
This reduces to: $\int y^{-2}dy=-\alpha \int x^{-2} dx=-y^{-1}=-\alpha x^{-1} + C$. This can be rewritten as:
$$y= \frac{1}{C+\frac{\alpha}{x}}, \text{for } x\neq 0$$
$$\text {undefined, for } x=0$$
The second part is because according to the solution, neither $x$, nor $y$ can be equal to zero, since this would mean division by zero. 
It seems we have to leave $y$ undefined for $x = 0$ because the intermediate solution contained a division by zero of both $y$ and $x$. Nevertheless, $Lim_{x\to0}(y)=0$, so this clearly suggests that we do not have to leave $y$ undefined at $x=0$, but can simply set $y(0)=0$. 
For example here is a graph of the function with $C=\alpha=1$:

My question is: How do we rigorously justify setting $y(0)=0$, even though in the solution we are dividing by both $y$ and $x$, so that it should be undefined at $y=0,x=0$?
 A: Your vector field is only well-defined for $x\neq 0$. The usual Cauchy Lipshitz condition is  not verified for $x=0$, so no guarantee a priori of existence/uniqueness of a solution through $x=0$.
It just happens so in this case, that the solutions from $x>0$ and $x<0$ may be joined as a $C^1$ function at $0$ when choosing the same const of integration on the two sides. For example if your ODE read $y'=-a y^2/x$, then the procedure wouldn't give rise to a $C^1$ function through 0, although it is continuous.
(there is btw a sign error in your solution).
A: I'm setting $\alpha=1$ since it doesn't affect this discussion.
$\dfrac{\mathrm{d}y}{\mathrm{d}x}=-\dfrac{y^{2}}{x^{2}}$ is undefined at $x=0$, but we may be able to get solutions with a point/removable discontinuity at $x=0$. One option is to say that with differential equations we want to fill in point discontinuities in our solutions; you probably wouldn't go wrong in most classes/applications with that mindset. I prefer to fix problems like these at the level of the original equation, saying “we should really be solving $x^{2}\dfrac{\mathrm{d}y}{\mathrm{d}x}=-y^{2}$”. 
Either way, by doing some standard steps under the assumptions that neither $y$ nor $x$ are zero, you can get a solution $y=-\dfrac{1}{C+\frac{1}{x}}$. Now, if $x^{2}\dfrac{\mathrm{d}y}{\mathrm{d}x}=-y^{2}$ then any solution must be zero at $x=0$. As you noted, ${\displaystyle \lim_{x\to0}}-\dfrac{1}{C+\frac{1}{x}}=0$, so this is a fine way to fill in the gap. 
Note that this doesn't require a piecewise definition: it just tells you that functions of the form $y=-\dfrac{x}{Cx+1}$ are solutions (even at $x=0$). 

This still leaves open the issue that we assumed $y\ne0$ to derive this family of solutions. Let's be more explicit about what we know. If a solution is nonzero somewhere, then since solutions should be continuous, it's nonzero on some tiny interval not containing $x=0$ (since $x=0$ implies $y=0$) around that point. On that tiny interval, our steps showed it follows a curve like $y=-\dfrac{x}{Cx+1}$, which it can follow until it hits a discontinuity or a place where $x$ or $y$ is zero (which only happens when both are). 
The only remaining case is if a solution is constantly $0$ (so it's not "nonzero somewhere"). Note that the constant $0$ function $y\equiv0$ does satisfy $x^{2}\dfrac{\mathrm{d}y}{\mathrm{d}x}=-y^{2}$, so that's another solution we missed. From one perspective, this is like saying $C=\infty$ is allowed in $y=-\dfrac{x}{Cx+1}$, but that sort of thinking won't necessarily get you a valid solution (or all of them) to every differential equation.
