Solve the differential equation $y'=y^2-x$ 
Solve the differential equation $$y'=y^2-x$$ with two different initial conditions: $y(0)= 1$ and $y(0)=0.5$.

My idea:
Suppose $y^2=t$ then $2yy'=t' \Rightarrow y'= \frac{t'}{2 \sqrt{t}}$
Given equation reduced by $\frac{t'}{2 \sqrt{t}}=t-x \Rightarrow t'=2 t\sqrt{t}-x \sqrt{t}$
How to solve from here  and how to use both initial conditions?
 A: With the substitution $$y=-\frac{u'}u,$$ the differential equation becomes linear:
$$u''=xu\tag1.$$ This equation is important in physics, that's why its solutions are well studied, the general solution is
$$u(x)=A\cdot\operatorname{Ai}(x)+B\cdot\operatorname{Bi}(x),$$ with special functions called Airy functions, and the constants can be determined from the initial condition
\begin{align}A\cdot\operatorname{Ai}(0)+B\cdot\operatorname{Bi}(0)&=1\\
A\cdot\operatorname{Ai}'(0)+B\cdot\operatorname{Bi}'(0)&=-y(0)\end{align}
or you can use the simpler equations
\begin{align}A&=\operatorname{Bi}'(0)+\operatorname{Bi}(0)\,y(0)\\B&=-\operatorname{Ai}'(0)-\operatorname{Ai}(0)\,y(0)\end{align}
(a constant factor of $u$ doesn't change our solution $y$). Since the functions $\operatorname{Ai}(x)$ and $\operatorname{Bi}(x)$ and their derivatives are implemented in Wolfram language, that's not more complicated than an expression using only standard functions. You can easily plot the solution corresponding to $y(0)=1$ in Wolfram Alpha and you'll see it has a singularity, a frequent phenomenon in non-linear differential equations. It's because of a zero of $u(x)$ at $x = 1.1073295812364012192025254$, and you can use Wolfram Alpha (or Mathematica) to determine that value, too.
Of course, you can also solve equation (1) by any numerical method. Since $u$ is an entire function, they'll all work very well. 
A: Solving the equation analytically is impossible by hand, as it involves a really complex expression involving the Bessel function. For such issues though, the solution approach via power series was implemented.
By the power series approach, assume that a solution will be proportional to : 
$$y(x) = \sum_{n = 0}^\infty a_n(x-x_0)^n$$
around an initial value point $x=x_0$.
For your specific case, since $x_0=0$ judging from the given IVPs, you will have to solve the following systems of equations for both cases : 
$$\begin{cases} y'=y^2-x \\ y(x) = \sum_{n = 0}^\infty a_nx^n \\ y(0) = 1 \end{cases}$$
$$\begin{cases} y'=y^2-x \\ y(x) = \sum_{n = 0}^\infty a_nx^n \\ y(0) = 1/2 \end{cases}$$
Doing so, you will have to calculate the derivative of the sum given to plug it in your differential equation and find out expressions and values for the constants $a_n$ :
$$y(x) = \sum_{n = 0}^\infty a_nx^n$$
$$y'(x) = \sum_{n = 0}^\infty na_nx^{n-1}$$
Can you proceed now ?
(Credits to Pedro as well for pointing out on the comments' section)
