# Using NDSolve for differential equations

I'm trying to numerically solve the following differential equations eqs for plotting the bifurcation diagram. However, under special values of the parameter u, I get a problem when calculating the differential equations with NDSolve. The corresponding codes are shown as follows

Initial definitions

\[Alpha] = 0.005;\[Beta] = 0.2;\[Gamma] = 10;kt = 100;\[Xi]b = 0.82;
\[Beta]L = {1.87510406871196117918., 4.69409113301900887678.,7.85475743823761498268., 10.99554073487734323228.,14.13716839104647178728., 17.27875953208826123958.,20.42035225210360409018.};
Subscript[\[Phi], m_][x] :=Cos[\[Beta]L[[m]]*x] - Cosh[\[Beta]L[[m]]*x] + ((Sin[\[Beta]L[[m]]] - Sinh[\[Beta]L[[m]]])/(Cos[\[Beta]L[[m]]] + Cosh[\[Beta]L[[m]]]))*(Sin[\[Beta]L[[m]]*x] - Sinh[\[Beta]L[[m]]*x]);
nig[fx_] :=N[Round[Chop[NIntegrate[fx, {x, 0, 1}, AccuracyGoal -> 10, WorkingPrecision -> 40, MaxRecursion -> 20], 10^-8], 10^-6], 6];


Differential equations

mm = 4;
eqs = Table[\!$$\*UnderoverscriptBox[\(\[Sum]$$, $$j = 1$$, $$mm$$]$$(\[Alpha]*nig[\(\*SubscriptBox[\(\[Phi]$$, $$i$$]\)[x]*D[$$\*SubscriptBox[\(\[Phi]$$, $$j$$]\)[x], {x, 4}]]*$$\*SubscriptBox[\(q$$, $$j$$]'\)[t])\)\)
+\!$$\*UnderoverscriptBox[\(\[Sum]$$, $$j = 1$$, $$mm$$]$$(nig[\(\*SubscriptBox[\(\[Phi]$$, $$i$$]\)[x]*D[$$\*SubscriptBox[\(\[Phi]$$, $$j$$]\)[x], {x, 4}]]*$$\*SubscriptBox[\(q$$, $$j$$]\)[t])\)\)
+\!$$\*UnderoverscriptBox[\(\[Sum]$$, $$j = 1$$, $$mm$$]$$(nig[\(\*SubscriptBox[\(\[Phi]$$, $$i$$]\)[x]*D[$$\*SubscriptBox[\(\[Phi]$$, $$j$$]\)[x], {x, 2}]]*\*SuperscriptBox[$$u$$, $$2$$]*$$\*SubscriptBox[\(q$$, $$j$$]\)[t])\)\)
-\!$$\*UnderoverscriptBox[\(\[Sum]$$, $$j = 1$$, $$mm$$]$$(nig[\(\*SubscriptBox[\(\[Phi]$$, $$i$$]\)[x]*D[$$\*SubscriptBox[\(\[Phi]$$, $$j$$]\)[x], {x, 2}]*\[Gamma]*$$(1 - x)$$]*$$\*SubscriptBox[\(q$$, $$j$$]\)[t])\)\)
+\!$$\*UnderoverscriptBox[\(\[Sum]$$, $$j = 1$$, $$mm$$]$$(nig[\(\*SubscriptBox[\(\[Phi]$$, $$i$$]\)[x]*D[$$\*SubscriptBox[\(\[Phi]$$, $$j$$]\)[x], x]]*2*\*SqrtBox[$$\[Beta]$$]*u*$$\*SubscriptBox[\(q$$, $$j$$]'\)[t])\)\)
+\!$$\*UnderoverscriptBox[\(\[Sum]$$, $$j = 1$$, $$mm$$]$$(nig[\(\*SubscriptBox[\(\[Phi]$$, $$i$$]\)[x]*D[$$\*SubscriptBox[\(\[Phi]$$, $$j$$]\)[x], x]]*\[Gamma]*$$\*SubscriptBox[\(q$$, $$j$$]\)[t])\)\)
+\!$$\*UnderoverscriptBox[\(\[Sum]$$, $$j = 1$$, $$mm$$]$$\*UnderoverscriptBox[\(\[Sum]$$, $$k = 1$$, $$mm$$]$$\*UnderoverscriptBox[\(\[Sum]$$, $$l = 1$$, $$mm$$]$$(\((\(\*SubscriptBox[\(\[Phi]$$, $$i$$]\)[x] /. x -> \[Xi]b)\)*$$(\(\*SubscriptBox[\(\[Phi]$$, $$j$$]\)[x] /. x -> \[Xi]b)\)*$$(\(\*SubscriptBox[\(\[Phi]$$, $$k$$]\)[x] /. x -> \[Xi]b)\)*$$(\(\*SubscriptBox[\(\[Phi]$$, $$l$$]\)[x] /. x -> \[Xi]b)\)*kt*$$\*SubscriptBox[\(q$$, $$j$$]\)[t]*$$\*SubscriptBox[\(q$$, $$k$$]\)[t]*$$\*SubscriptBox[\(q$$, $$l$$]\)[t])\)\)\)\)
+\!$$\*UnderoverscriptBox[\(\[Sum]$$, $$j = 1$$, $$mm$$]$$(nig[\(\*SubscriptBox[\(\[Phi]$$, $$i$$]\)[x]*$$\*SubscriptBox[\(\[Phi]$$, $$j$$]\)[x]]*$$\*SubscriptBox[\(q$$, $$j$$]''\)[t])\)\) == 0, {i, mm}];
EQs = Chop[Expand[Join[eqs]], 10^-8]


Initial conditions

{q101, q102, q201, q202, q301, q302, q401, q402} = {0.001, 0, 0, 0, 0,
0, 0, 0};
ICs = {Subscript[q, 1] == q101, Subscript[q, 1]' == q102, Subscript[q, 2] == q201, Subscript[q, 2]' == q202, Subscript[q, 3] == q301, Subscript[q, 3]' == q302, Subscript[q, 4] == q401, Subscript[q, 4]' == q402};

u = 8.8;
s0 = NDSolve[Join[EQs, ICs], {Subscript[q, 1], Subscript[q, 2], Subscript[q, 3], Subscript[q, 4]}, {t, 0, 1000}, Method -> "ImplicitRungeKutta"]


The problem is that when the parameter value u is lower than 8.8, like u=8.7 or even u=6.5, the above-mentioned codes are efficient for the calculation in the interval $$u\in\left(6,8.7\right)$$. But when the parameter u beyonds 8.7, such as u=8.8, 8.9, and 9, etc, I get the following warning:

NDSolve::mxst: Maximum number of 216267 steps reached at the point t == 83.2334874133867

Concerning this warning, I use

• MaxSteps -> 1000000
• WorkingPrecision -> 40

But Mathematica still gives me the warning. When selecting MaxSteps -> Infinity, Mathematica will keeps running for a long time. And then, I stops it.

## Question

I would like to numerically solve the differential equations eqs with $$u\in\left(6,9.5\right)$$. So I am wondering:

• Is there a specific method for numerically solving the differential equations? It may be not the "ImplicitRungeKutta" or even the "NDSolve".
• I succeeded in calculating these equations when choosing mm=2, where mm is the number of 2-order differential equations considered. I suspect that the appearance of the warning NDSolve::mxst: Maximum number of $$\times\times$$ steps reached at the point t == $$\times\times$$ may have something to do with the value of mm, but I am not sure. Is there any other way I can do that when mm>4?

Thanks for help and insights in advance!

## Supplement

I have noticed MarcoB's suggestions. Some related replies could be seen in the comment section. Here, I add some changes to the subscript format of the code.

Initial definitions

\[Phi][m_, x_] := Cos[\[Beta]L[[m]]*x] - Cosh[\[Beta]L[[m]]*x] + ((Sin[\[Beta]L[[m]]] - Sinh[\[Beta]L[[m]]])/(Cos[\[Beta]L[[m]]] + Cosh[\[Beta]L[[m]]]))*(Sin[\[Beta]L[[m]]*x] - Sinh[\[Beta]L[[m]]*x]);


Differential equations

eqs = Table[\!$$\*UnderoverscriptBox[\(\[Sum]$$, $$j = 1$$, $$mm$$]$$(\ [Alpha]*nig[\[Phi][i, x]*D[\[Phi][j, x], {x, 4}]]*D[q[j, t], t])$$\)
+ \!$$\*UnderoverscriptBox[\(\[Sum]$$, $$j = 1$$, $$mm$$]$$(nig[\[Phi][i, x]*D[\[Phi][j, x], {x, 4}]]*q[j, t])$$\)
+ \!$$\*UnderoverscriptBox[\(\[Sum]$$, $$j = 1$$, $$mm$$]$$(nig[\[Phi][i, x]*D[\[Phi][j, x], {x, 2}]]*\*SuperscriptBox[\(u$$, $$2$$]*q[j, t])\)\)
- \!$$\*UnderoverscriptBox[\(\[Sum]$$, $$j = 1$$, $$mm$$]$$(nig[\[Phi][i, x]*D[\[Phi][j, x], {x, 2}]*\[Gamma]*\((1 - x)$$]*q[j, t])\)\)
+ \!$$\*UnderoverscriptBox[\(\[Sum]$$, $$j = 1$$, $$mm$$]$$(nig[\[Phi][i, x]*D[\[Phi][j, x], x]]*2*\*SqrtBox[\(\[Beta]$$]*u*D[q[j, t], t])\)\)
+ \!$$\*UnderoverscriptBox[\(\[Sum]$$, $$j = 1$$, $$mm$$]$$(nig[\[Phi][i, x]*D[\[Phi][j, x], x]]*\[Gamma]*q[j, t])$$\)
+ \!$$\*UnderoverscriptBox[\(\[Sum]$$, $$j = 1$$, $$mm$$]\
(\*UnderoverscriptBox[$$\[Sum]$$, $$k = 1$$, $$mm$$]$$\*UnderoverscriptBox[\ (\[Sum]$$, $$l = 1$$, $$mm$$]$$(\((\[Phi][i, x] /. x -> \[Xi]b)$$*$$(\[Phi] [j, x] /. x -> \[Xi]b)$$*$$(\[Phi][k, x] /. x -> \[Xi]b)$$*$$(\[Phi][l, x] /. x -> \[Xi]b)$$*kt*q[j, t]*q[k, t]*q[l, t])\)\)\)\)
+ \!$$\*UnderoverscriptBox[\(\[Sum]$$, $$j = 1$$, $$mm$$]$$(nig[\[Phi][i, x]*\[Phi][j, x]]*D[q[j, t], {t, 2}])$$\) == 0, {i, mm}];

EQs = EQs /. {q[1, t] -> q1[t], q[2, t] -> q2[t], q[3, t] -> q3[t], q[4, t] -> q4[t], \!$$\*SuperscriptBox[\(q$$, TagBox[RowBox[{"(", RowBox[{"0", ",", "1"}], ")"}],Derivative],MultilineFunction->None]\)[1, t] -> q1'[t], \!(\*SuperscriptBox[$$q$$, TagBox[RowBox[{"(", RowBox[{"0", ",", "1"}], ")"}],Derivative],MultilineFunction->None]\)[2, t] -> q2'[t], \!$$\*SuperscriptBox[\(q$$, TagBox[RowBox[{"(", RowBox[{"0", ",", "1"}], ")"}],Derivative],MultilineFunction->None]\)[3, t] -> q3'[t], \!$$\*SuperscriptBox[\(q$$, TagBox[RowBox[{"(", RowBox[{"0", ",", "1"}], ")"}],Derivative],MultilineFunction->None]\)[4, t] -> q4'[t], \!$$\*SuperscriptBox[\(q$$, TagBox[RowBox[{"(", RowBox[{"0", ",", "2"}], ")"}],Derivative],MultilineFunction->None]\)[1, t] -> q1''[t], \!$$\*SuperscriptBox[\(q$$, TagBox[RowBox[{"(", RowBox[{"0", ",", "2"}], ")"}],Derivative],MultilineFunction->None]\)[2, t] -> q2''[t], \!$$\*SuperscriptBox[\(q$$, TagBox[RowBox[{"(", RowBox[{"0", ",", "2"}], ")"}],Derivative],MultilineFunction->None]\)[3, t] -> q3''[t], \!$$\*SuperscriptBox[\(q$$, TagBox[RowBox[{"(", RowBox[{"0", ",", "2"}], ")"}],Derivative],MultilineFunction->None]\)[4, t] -> q4''[t]}


Initial conditions

q1 == 0.001, q1' == 0, q2 == 0, q2' == 0, q3 == 0, q3' == 0, q4 == 0, q4' == 0

• Did you notice that your initial conditions are defined as ICs, but then in your NDSolve call you use an undefined ICS1? Is that a typo? Also, your betaL values are defined with 8 digits of precision; it does not make sense to ask for a WorkingPrecision of 40 later on. Have you tried removing all the WorkingPrecision, and explicit Method calls? What happens if you let NDSolve choose the method automatically? Also, Subscripts are mostly a formatting convenience, but they are best avoided in calculation, where they might give rise to sneaky problems. – MarcoB Dec 23 '19 at 16:08
• @MarcoB thanks. Firstly, I would change the ICS1 to ICS later. It was right in my Mathematica codes before. Secondly, when the "WorkingPrecision->40" is neglected, the warning "the numerical integral converges too slowly" will happends. Changing the Subscripts and removing the WorkingPrecision or the explicit method, the warning "NDSolve::mxst: Maximum number of ×××× steps reached at the point t == ××" still has. Could it be something else？ – kuzb Dec 24 '19 at 2:50
• @MarcoB I have tried to numerically solve the differential equations with mm=2 using the ode45 of Matlab. At u=8.8, 8.9, and 9, Matlab could easily solve it. The possibility of the problem may comes from the Method. ode45-NDSolve – kuzb Dec 24 '19 at 13:59