The equation $$ y' = 2\sqrt{y} $$ is a first order, nonlinear, autonomous, constant coefficient ordinary differential equation (DE). Since it is first order and autonomous, it is susceptible to separation of variables, and WA says as much. So we obtain
$$ 2 \sqrt{y} + C_1 = \int \frac{\mathrm{d}y}{\sqrt{y}} = \int 2 \,\mathrm{d}x = 2x + C_2 $$
Assuming $y \geq 0$ for all $x$, this can be simplified to $y = x^2 + 2Cx + C^2, y \geq 0$. Since WA (and Mathematica) ignore constants of integration in antiderivatives, WA only writes "$y = x^2$". However, this hides an important property of autonomous systems: horizontally translating a solution produces another solution. We can exhibit that translatability in our solution:
$$ y = x^2 + 2Cx + C^2 = (x+C)^2 $$
and plugging back into the original equation,
$$ 2(x+C) = y' = 2 \sqrt{y} = 2\sqrt{(x+C)^2} = 2|x+C| \text{,} $$
which requires $x \geq -C$ (or, equivalently, $y' \geq 0$). So separation of variables has given us the family of solutions $y = (x+C)^2, x \geq -C$.
The original DE is nonlinear so there is no promise that finding a single family of solutions exhausts the solution space. For any homogeneous linear DE, $y = 0$ is a solution. This system is not homogeneous and linear, but it is autonomous and the right-hand side, $2\sqrt{y}$, has a zero at $y = 0$, so we are lucky enough to also have the additional solution $y = 0$.
From what we have done so far, there could well be more isolated solutions or solution families we have not found. Nonlinear DEs do not come with the existence and uniqueness guarantees that linear DEs do and we would need some sort of uniqueness to be certain we had exhausted the set of solutions.
However, what we have done shows the following: you can glue the $y =0$ solution on $(-\infty, -C)$ to the half-parabola $(x+C)^2, x \geq -C$ and get another solution. That is, the solution $y = 0$ can, at any time, suddenly jump off the $x$-axis along a positive half-parabola, all the while satisfying the given DE. If that seems impossible, note that the equation is satisfied by $y' = y = 0$ and both the $y =0$ solution and the half-parabola have $y' = y = 0$ at the gluing point. We might think that there should be some difference in higher derivatives. There is, but that doesn't change any of the analysis so far.
One result is that the initial condition $y(0) =0$ does not pick out a single solution. Uncountably many solutions pass through the origin, differing only in the point on the nonnegative $x$-axis that is the
vertex of their parabola, except for the one additional solution that never leaves the $x$-axis. This is very different from linear DEs, where only one solution passes through a particular initial condition.
If we differentiate the given equation, we obtain
$$y'' = y' y^{-1/2} \text{.}$$
The right-hand side is undefined at $y = 0$, so we cannot trust that this equation tells us anything about solutions to the original equation when $y = 0$. If we blindly substitute the original equation, we obtain $y'' = \frac{2\sqrt{y}}{\sqrt{y}} = 2, y \neq 0$. Of course, this is true for all the half-parabola solutions. This bad behaviour at $y = 0$ tells us something odd/interesting should be occurring at/near the solution $y = 0$ and that this oddness/interestingness is caused by the square root in $\sqrt{y}$.
Let's take a brief detour to discuss WolframAlpha's and Mathematica's handling of branched functions. These pieces of software are variously oblivious, silently switching branches (the WA solution you reference), or meticulously careful about branches. There are two branches of the (real) square root. The nonnegative branch is denoted "$\sqrt{\ }$", the square root function, and the other branch is obtained by negation. Unfortunately, these pieces of software may ignore our notation and/or switch branches without warning. If we study the DE, but switch the branch of the square root, we analyze $$y' = -2\sqrt{y} \text{.}$$
Doing so, we still get $y = 0$ and we get the left-halves of unit parabolas with vertices on the $x$-axis -- the mirror images of the families of parabolas solving the first equation. So, for this choice of square root, the general solutions are descending half-parabolas that turn into the $y= 0$ solution when they reach their vertices on the $x$-axis. Perhaps this seems less weird than the "spontaneous jumping away" from the $x$-axis that we got with the original equation, but both behaviours are and should seem equally weird.
But this highlights how WA's whole parabola solution has an interpretation:
- pick a positively opening unit parabola with vertex on the $x$-axis,
- on the left half of the parabola, pick the nonpositive branch of the square root and see that the parabola satisfies the DE,
- at the vertex of the parabola, observe that the equation is satisfied with either choice of the branch of the square root, and
- on the right half of the parabola, pick the nonnegative branch of the square root and see that the parabola continues to satisfy the DE.
This is, of course, bogus, because the notation "$\sqrt{y}$" denotes the square root function (which is the nonnegative branch), not "some freely changing branch of the square root". That is, once we have all agreed that "$\sqrt{y}$" means the nonnegative branch of the square root, then WA's $y=x^2$ is wrong. "In reality", only the ascending halves of these parabolas are solutions to the given DE and these solutions are obtained for "generic" initial conditions. Only special initial conditions give the $y=0$ solution.
This discussion of branches and $y''$ has ranged a bit, but it has been preparation for answering your second question.
The $y^{-1/2}$ in the second derivative, $y'' = y' y^{-1/2}$, has very rapidly changing behaviour for small, positive $y$, but the parabola solutions are very slowly varying there. This is a (not entirely satisfactory) signal that the DE can be stiff for certain numerical methods. It is worth pointing out that ode45 is used to "Solve nonstiff differential equations", so we should generically expect that this solver has trouble with stiffness. There are aspects of stiffness present in the transition from $y = 0$ to $y > 0$ so we should expect that this solver may be unable to reproduce that transition. (In fact, I don't expect any numerical solver to switch from the $y=0$ solution to a half-parabola solution unless it is so numerically noisy that it can't follow the $y=0$ solution, which would make it largely useless for any other problem.)
What happens "in reality" is that ode45 happily follows half-parabola solutions when given an initial condition off the $x$-axis. That means that for almost all initial conditions (i.e., $y(x_0) = y_0 > 0$ for any $x_0 \in \Bbb{R}$), ode45 happily follows a half-parabola solution. But for the special initial conditions on the boundary of these (i.e., $y(x_0) = 0$ for any $x_0 \in \Bbb{R}$), ode45 follows the $y = 0$ solution.
So, if anything, the $y=0$ solution is the weird solution that traps the solver. ode45 should follow the right half of a WA solution if given an initial condition that doesn't trap it in the $y=0$ solution.
I'm somewhat curious what ode45 does with the "other square root branch" equation when given various initial conditions above the $x$-axis. I expect it follows their parabolas down to the $x$-axis but I'm not sure if it's ever lucky enough to subsequently follow the $y = 0$ solution or is routinely unlucky enough to overshoot to $y < 0$ and abort due to non-real evaluation of $\sqrt{y}$.
A brief comment on nonnegativity of $y$ that can be used to avoid the non-real evaluation problem mentioned above. One way to enforce that $\sqrt{y}$ always exists is to replace $y$ with a function that is always nonnegative. One way is to write $y = u^2$, where $u$ is a new real-valued dependent variable. The resulting DE is $(u^2)'=2|u|$. Passing $(u^2)'=2u$ to WA, the solutions are $u = 0$ and $u = x+C$, from which $y = 0$ and $y = (x+C)^2$. Passing $(u^2)'=-2u$ to WA, the solutions are $u = 0$ and $y = -x + C$, from which $y = 0$ and $y=(x+C)^2$. So WA can give the exceptional $y = 0$ solutions as well as misrepresenting the family of half-parabola solutions. We have to be smart enough to figure out how to glue the solutions to the two equations together so that the resulting solutions satisfy the version with absolute values -- when $u \geq 0$, $2|u| = 2u$ and when $u \leq 0$, $2|u| = -2u$, effectively expressing the choice of branch of the square root as choice of sign of $u$.
There are other ways to enforce nonnegativity of $y$. For instance, set $y = \mathrm{e}^v$ and remember that $y(x) = \mathrm{e}^{v(x) + 2 \pi \mathrm{i} k}$ for any integer $k$. The DE becomes $$\frac{\mathrm{d}}{\mathrm{d}x} \mathrm{e}^{v(x) + 2 \pi \mathrm{i} k} = 2 \mathrm{e}^{\frac{1}{2}v(x) + \pi \mathrm{i} k}, k \in \Bbb{Z} \text{.} $$ (Here we encode the choice of branch of the square root in the expression "$\pi \mathrm{i} k$" by our choice of even or odd $k$.)
The difficulty with this method is the initial condition, $v(0) = -\infty$, which we might recognize is ... problematic for a numerical solver. Even if a solver could somehow start at $-\infty$, it can only subsequently take finite steps, arriving at $-\infty$. So we have, perhaps, found another way to express that a $y(x_0) = 0$ initial condition traps solvers on the solution $y = 0$, but any positive initial condition, $y(x_0) = y_0 > 0$ yields a completely adequate numerical solution along a half-parabola.
There is no one-size-fits-all way to handle nonlinear DEs. So there is no generic answer to your question about recognizing "bad" solutions. When facing a new, unfamiliar nonlinear DE, expect that some solutions may be found by applying techniques you know, some solutions may require techniques you don't know, and you have no way to know if you've found all the solutions. The longer you work with nonlinear DEs and watch/read what others do with them, the more "techniques you know" you will have for addressing a new one. You may find that certain nonlinear equations do have complete solution techniques, e.g., Bernoulli differential equations, but you will learn that these are exceptional -- generic nonlinear DEs are complicated.
