# How to solve second order nonlinear differential equation?

sol = NDSolve[{y''[x] + (y[x])^2 == 1, y[0] == 2, y'[0] == 3}, y[x], {x, 0, 30}]
g = Plot[Evaluate[{y[x], y'[x], y''[x]} /. sol], {x, 0, 30}, PlotStyle -> Automatic]


I'm trying to coding this equation $$y''+ y^2 =1$$ with initial value $$y(0)= 2$$, and $$y'(0)= 3$$ using NDSolve. Anyone could tell me what is the problem above in the code? Thank you.

• It does solve the system, but I get a message which says that there is a singularity at x = 3.95. Oct 12, 2018 at 5:13
• Also, by changing the initial condition to y'[0] == 0 gives a nice periodic solution. Is this what you want ? Oct 12, 2018 at 5:17
• I need it to give me a graph, but I tried old method I learned before, and it didn't work. Oct 12, 2018 at 5:22
• How to change the initial value for y ' [0]=3 to y ' [0]=0 when the given initial in the problem is this one "y ' [0]=3" Oct 12, 2018 at 5:25
• I am not saying that we should change it. I am just pointing out that with y'[0] = 3 we get stiffness issues, which goes away with y'[0] = 0. Understanding of the physics of the problem may explain this behaviour. Oct 12, 2018 at 5:33

The desired solution can be obtained semi-analytically.

s = DSolveValue[y''[x] + (y[x])^2 == 1, y[x], x] /. {C[1] -> c1, C[2] -> c2}
(* -(-1)^(2/3) 6^(1/3)
WeierstrassP[(-(1/6))^(1/3) (c1 + x), {2 (-1)^(2/3) 6^(1/3), c2}] *)


Next, determine the constants by applying the boundary conditions.

bc1 = s /. x -> 0;
bc2 = D[s, x] /. x -> 0;
sc = FindRoot[{bc1 == 2, bc2 == 3}, {c1, I}, {c2, 1}] // Chop
(* {c1 -> -0.523953 + 1.79242 I, c2 -> -10.3333} *)


and plot.

Plot[Chop[s /. sc], {x, 0, 10}, ImageSize -> Large]


The solution is periodic in x. In light of the singularities in the solution, it is not surprising that NDSolve failed for x near 4.

The period can be increased by decreasing y'[0] and becomes infinite for y'[0] == 0. A WorkingPrecision of at least 45 is needed in that case.

May not be exactly what OP wants but still:

sol = NDSolve[{y''[x] + (y[x])^2 == 1, y[0] == 2, y'[0] == 0},
y, {x, 0, 30}]

g = Plot[Evaluate[{y[x], y'[x], y''[x]} /. sol], {x, 0, 30},
PlotStyle -> Automatic]


Gives the graphs as required. With y'[0] == 3 as given by OP there seems to be a singularity/stiffness at x = 3.95

just a supplement to @Lotus answer:

Using ParametricNDSolveValue

Y = ParametricNDSolveValue[{y''[x] + (y[x])^2 == 1, y[0] == 2,y'[0] == ys}, y, {x, 0, 30}, ys]


you could do some parameter studies varying ys=y'[0]

Show[{
Plot[Table[Y[p][x], {p, .1, 3. , .1}] , {x, 0, 25},PlotRange -> {-5, 5}] ,
Plot[ Y[0][x] , {x, 0, 25}, PlotStyle -> {Thickness[.01], Red}],
Plot[ Y[3][x] , {x, 0, 25}, PlotStyle -> { Green },PlotRange -> {-5, 5}]
}]


which show, that only y'[0]==0 gives convergent solution! (The green solution is y'[0]==3 OP asked for )

• Increasing WorkingPrecision will give a different plot for y'[0]==0. Oct 12, 2018 at 15:34
• @bbgodfrey Thanks, trying ParametricNDSolveValue[... , WorkingPrecision -> 50] Plot[Y[3][x], {x, 0, 10}, PlotStyle -> {Green}, PlotRange -> {-5, 5} ] doesn't finish ... Oct 12, 2018 at 15:54
• Works fine for me with NDSolveValue and Mathematica 11.3. Solution is flat at about -1 for large x, which agrees with the symbolic solution. Oct 12, 2018 at 16:05
• @bbgodfrey Thanks, MMA 11.0.1. NDSolve gives the solution in the range 0<x<3.95 Oct 12, 2018 at 19:25