Using NDSolve for differential equations

Question

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.8751040687119611791`8., 4.6940911330190088767`8.,7.8547574382376149826`8., 10.9955407348773432322`8.,14.1371683910464717872`8., 17.2787595320882612395`8.,20.4203522521036040901`8.};
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][0] == q101, Subscript[q, 1]'[0] == q102, Subscript[q, 2][0] == q201, Subscript[q, 2]'[0] == q202, Subscript[q, 3][0] == q301, Subscript[q, 3]'[0] == q302, Subscript[q, 4][0] == q401, Subscript[q, 4]'[0] == 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.