Solve the differential equation $y'=y^2-x$

Question

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?

Answer 1

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)

Answer 2

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.