# How to NDSolve this kind of ODE?

For an ODE like this:$(1-y)y'+y^2=0$ with the initial condition $y(1)=1$, how to solve it numerically? I know this equation can be solved analytically by DSolve. In fact, my equation is more complicated than this, I have to solve it numerically. Using NDSolve directly,

NDSolve[{(1 - y[x])*y'[x] + y[x]^2 == 0, y[1] == 1}, y, {x, 1, 5}]


it will display error messages:

Power::infy: Infinite expression 1/0. encountered.
NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 1..


I guess this problem happens because the initial condition just makes the coefficient of y'[x] be zero. So my question is how to overcome this problem?

• The derivative diverges at the point $x=1$, and it's hard to deal with divergencies numerically. Perturbing the initial condition slightly removes the problem, e.g. NDSolve[{(1 - y[x])*y'[x] + y[x]^2 == 0, y[1] == 1.01}, y, {x, 1, 5}] – Marius Ladegård Meyer Nov 22 '16 at 9:57

One can here introduce another dependent variable: z[x]->y[x] - y[x]^2/2and express your equation in terms of this variable:

    ss = NDSolve[{z'[x] + (1 - Sqrt[1 - 2 z[x]])^2 == 0, z[1] == 1/2},
z[x], {x, 1, 10}][[1, 1]]
(*  z[x] -> InterpolatingFunction[{{1., 10.}}, <>][x]  *)


which is nicely solved:

Plot[1 - Sqrt[1 - 2 z[x]] /. ss, {x, 1, 10},
AxesLabel -> {Style["x", 18, Italic], Style["y", 18, Italic]}]


yielding

Have fun!

• This transformation actually will lead to 2 equations: ((1 - y[x])*y'[x] + y[x]^2 == 0 /. Solve[z[x] == y[x] - y[x]^2/2, y[x]] /. (y[x] -> a_) :> (y -> (Function[x, #] &@a))), which correspond to 2 solutions of the original equation :) – xzczd Nov 22 '16 at 11:12
• @ xzczd of course, but since this is evident, the OP can figure it out himself and decide, what solution is he interested in. – Alexei Boulbitch Nov 22 '16 at 12:26

Your specific example is actually special, it has 2 solutions, and DSolve can only find one of them. To find both of the solutions, we can modify the equation from a equation of $y(x)$ to a equation of $x(y)$:

$$\frac{1-y}{x'(y)}+y^2=0$$

Then DSolve and NDSolve can both handle the equation without difficulty:

asolinverse = x /. First@DSolve[{(1 - y)*1/x'[y] + y^2 == 0, x[1] == 1}, x, y]
(* Function[{y}, (1 + y Log[y])/y] *)
nsolinverse =
x /. First@NDSolve[{(1 - y)*1/x'[y] + y^2 == 0, x[1] == 1}, x, {y, 10^-3, 100}]

ParametricPlot[{nsolinverse[y], y}, {y, 10^-3, 10}, AspectRatio -> 1/GoldenRatio]
`

• your method is also very good. But I think Boulbitch's method is more universal. For more complicated ODEs, it is not easy to modify the equation from an equation of y(x) to an equation of x(y) . – Mark_Phys Nov 23 '16 at 10:46