# Any ideas on how to use NDSolve to solve this nonlinear 2nd order ODE

I can't get to work this code to solve a nonlinear second order boundary value problem.

ClearAll["Global*"];
a = 1;
b = 2;
q1[t_] := {Sin[2*Pi*t], Cos[2*Pi*t]}
q2[t_] := {a*Sin[2*Pi*t], b*Cos[2*Pi*t]}
f[x_, y_] := x . y
eqn = 0.5*Derivative[2][y][t]*f[q2[y[t]], q1[t]] -
3*Derivative[1][y][t]^2*
f[Derivative[1][q2][y[t]], q1[t]] +
Derivative[1][y][t]*f[q2[y[t]], Derivative[1][q1][t]] +

4*Derivative[1][y][t]^2*Sqrt[Derivative[1][y][t]]*
f[q2[y[t]], q2[y[t]]] == 0;
bc = {y[0] == 0, y[1] == 1};
sol = NDSolve[{eqn, bc}, y, {t, 0, 1},
Method -> "StiffnessSwitching"]
Plot[Evaluate[y[x] /. sol], {x, 0, 1}, PlotRange -> All]


I get immediately the error message

(1) "NDSolve::ndsz: At t == 0.2499999999990644, step size is effectively zero; singularity or stiff system suspected."

When I add "AccuracyGoal -> 18, PrecisionGoal -> 18" to this code, I get the message

(2) "NDSolve::mxst: Maximum number of 79973 steps reached at the point t == 0.2499927591408232."

When I reduce to "AccuracyGoal -> 8, PrecisionGoal -> 8", I get message (1).

One difficulty in solving this ODE may be the square root in the equation. So, I converted it into a nonlinear system of first order equations without the square root using the transformation $$u=\sqrt{y'}$$. Then, $$u^2 = y'$$ and $$2uu'=y''$$. Also, let $$v=y$$, then $$v'=u^2$$. The 1st order system is: $$\{u'-u^3f_1(t,v)+uf_2(t,v)+u^4f_3(t,v)=0,v'=u^2\}$$, where $$f_1(t,v)=q_2'(y).q_2(t)/q_2(y).q_1(t)$$, $$f_2(t,v)=q_2(y).q_1'(t)/q_2(y).q_1(t)$$, $$f_3(t,v)=q_2'(y).q_2(t)/q_2(y).q_1(t)$$

The code is:

s = NDSolve[{Derivative[1][u][t] -
u[t]^3*(3*(f[Derivative[1][q2][v[t]], q1[t]]/
f[q2[v[t]], q1[t]])) +
u[t]*(f[q2[v[t]], Derivative[1][q1][t]]/
f[q2[v[t]], q1[t]]) + u[t]^4*
(4*(f[Derivative[1][q2][v[t]], q2[v[t]]]/
f[q2[v[t]], q1[t]])) == 0,
Derivative[1][v][t] - u[t]^2 == 0, v[0] == 0, v[1] == 1},
{u, v}, {t, 0, 1}, Method -> "StiffnessSwitching",
AccuracyGoal -> 18, PrecisionGoal -> 18]
Plot[Evaluate[{v[t], u[t]}], {t, 0, 1}, PlotRange -> All]


I get the message, "The scaled boundary value residual error of 4.714045207910324$$$$*^7 indicates that the boundary values are not satisfied to specified tolerances. Returning the best solution found." It returns constant solutions.

Many nonlinear ODEs have no solution and this is such a case, as can be seen from

ff = ParametricNDSolveValue[{eqn, y[0] == 0, y'[0] == y0}, y[1], {t, 0, 1}, {y0}];
Plot[ff[w], {w, 0, 4}, AxesLabel -> {y'[0], y[1]},
LabelStyle -> {12, Bold, Black}, PlotRange -> All]


For no value of y'[0] does y[1] equal 1. Here are some typical solutions.

gg = ParametricNDSolveValue[{eqn, y[0] == 0, y'[0] == y0}, y[t],
{t, 0, 1}, {y0}];
Plot[Evaluate@Table[gg[i], {i, {.01, .1, 1, 4}}], {t, 0, 1},
AxesLabel -> {t, y[t]}, LabelStyle -> {12, Bold, Black}, PlotRange -> All]


• Thanks very much. These solutions look very promising. I can't understand, though, why you can't reach y(1)=1. Do you? Commented May 30 at 13:33
• The behavior of nonlinear ODEs is difficult to predict, much less to explain. Clearly, the solution has a strong contractor, because all values of y'[0] soon converge to approximately the same solution. And that solution does not reach 1. Commented May 30 at 20:53