nonlinear equation numerical solve method for stable calculation

I clearify my question. I write a code that diverges strangly after 8.4 second. but my system is completely stable and I think Mathematica's solution method is not true. Please help me for find complete and correct result. If i choose a time more than diverge time my code shows an error!!!

NDSolve::mconly: For the method IDA, only machine real code is available. Unable to continue with complex values or beyond floating-point exceptions.

Or using some of initial condition can make problem in solve my equations!!! My system is completely stable for every initial condition but result of my code is not true and stable for some inputs and initial conditions. This code is solution of a set of nonlinear differential equation that shows dynamical behavior of a stable mechanical system. Complexity of my system is more than this equations. I need a solving method for a very complicated dynamics behavior equations. Is Mathematica a proper software for my work? Or i should use other software such as maple or Matlab? in your opinion can I use this software Predetermined functions like ndsolve(Mathematica), ODE45(Matlab), ... or i should write solving codes completely by my own.

ClearAll[t, L, L1, L2, L3, θ1, θ2, θ3,
θ4, θ5, θ6, eq1, eq2, eq3, eqs]
g = 9.81;
R = 0.3;
h = 0.1;

thickness = 0.01;

L1 = 0.5;
L2 = 0.6;
L3H = 0.1;
L3V = 0.05;
L4 = 0.5;
L5 = 0.4;
L6V = L3V;
L6H = L3H;
e1 = L6H/4;
e2 = L6H/2;

LC1 = L1/2;
LC2 = L2/2;
LC3V = L3V/2;
LC3H = L3H/2;
L3H2 = L3H*2/3;
LC4 = L4/2;
LC5 = L5/2;
LC6H = L6H/2;
LC6V = L6V/2;

(*Mass*)
m1 = 1;
m2 = 1;
m3 = 1;
m4 = 1;
m5 = 1;
m6 = 1;

(*Moment of inertia*)
I1 = 1;
I2 = 1;
I3 = 1;
I4 = 1;
I5 = 1;
I6 = 1;

(*Viscos damping factors*)
C1 = 2;
C2 = 2;
C3 = 2;
C4 = 2;
C5 = 2;
C6 = 2;

RDF = 1/2*C1*(θ1'[t])^2 + 1/2*C2*(θ2'[t])^2 +
1/2*C3*(θ3'[t])^2 + 1/2*C4*(θ4'[t])^2 +
1/2*C5*(θ5'[t])^2 + 1/2*C6*(θ6'[t])^2;

(*Input torques*)
T1 = 0;
T2 = 0;
T3 = 0;
T4 = 5 Sin[t];
T5 = 5 - 5*Sin[t];
T6 = 5 Sin[t];

xe[t_] := -L4*Sin[θ4[t]] -
L5*Cos[Pi/2 - θ4[t] + θ5[t]] -
L6V*Cos[Pi/2 - θ4[t] + θ5[t] + θ6[t]] -
L6H*Sin[Pi/2 - θ4[t] + θ5[t] + θ6[t]] +
L3V*Cos[Pi/2 - θ4[t] + θ5[t] + θ6[t]];
ye[t_] := -(R - h) - L4*Cos[θ4[t]] -
L5*Sin[Pi/2 - θ4[t] + θ5[t]] -
L6V*Sin[Pi/2 - θ4[t] + θ5[t] + θ6[t]] +
L6H*Cos[Pi/2 - θ4[t] + θ5[t] + θ6[t]] +
L3V*Sin[Pi/2 - θ4[t] + θ5[t] + θ6[t]];

a = L1;
b = R;
c[t_] := (xe[t]^2 + ye[t]^2 - R^2 - L1^2 - L2^2)/(2 L2);

θ2[t_] := 2 ArcTan[(b + Sqrt[b^2 + a^2 - c[t]^2])/(a + c[t])]
θ1[t_] :=
ArcTan[(-(L1 + L2*Cos[θ2[t]])*
xe[t] - (R + L2*Sin[θ2[t]])*
ye[t]), (+(R + L2*Sin[θ2[t]])*
xe[t] - (L1 + L2*Cos[θ2[t]])*ye[t])]
θ3[t_] := θ6[t] + θ5[t] +
Pi/2 - θ4[t] - (θ1[t] + θ2[t]);

(*Exoskeleton*)
XG1 = +R*Sin[θ1[t]] - LC1*Cos[θ1[t]];
YG1 = -R*Cos[θ1[t]] - LC1*Sin[θ1[t]];
XG2 = +R*Sin[θ1[t]] - L1*Cos[θ1[t]] -
LC2*Cos[θ1[t] + θ2[t]];
YG2 = -R*Cos[θ1[t]] - L1*Sin[θ1[t]] -
LC2*Sin[θ1[t] + θ2[t]];
XG3 = +R*Sin[θ1[t]] - L1*Cos[θ1[t]] -
L2*Cos[θ1[t] + θ2[t]] -
LC3V*Cos[θ1[t] + θ2[t] + θ3[t]] +
LC3H*Sin[θ1[t] + θ2[t] + θ3[t]];
YG3 = -R*Cos[θ1[t]] - L1*Sin[θ1[t]] -
L2*Sin[θ1[t] + θ2[t]] -
LC3V*Sin[θ1[t] + θ2[t] + θ3[t]] -
LC3H*Cos[θ1[t] + θ2[t] + θ3[t]];

XG4 = -LC4*Sin[θ4[t]];
YG4 = -(R - h) - LC4*Cos[θ4[t]];
XG5 = -L4*Sin[θ4[t]] -
LC5*Cos[Pi/2 - θ4[t] + θ5[t]];
YG5 = -(R - h) - L4*Cos[θ4[t]] -
LC5*Sin[Pi/2 - θ4[t] + θ5[t]];
XG6 = -L4*Sin[θ4[t]] -
L5*Cos[Pi/2 - θ4[t] + θ5[t]] -
LC6V*Cos[Pi/2 - θ4[t] + θ5[t] + θ6[t]] -
LC6H*Sin[Pi/2 - θ4[t] + θ5[t] + θ6[t]];
YG6 = -(R - h) - L4*Cos[θ4[t]] -
L5*Sin[Pi/2 - θ4[t] + θ5[t]] -
LC6V*Sin[Pi/2 - θ4[t] + θ5[t] + θ6[t]] +
LC6H*Cos[Pi/2 - θ4[t] + θ5[t] + θ6[t]];

ω1 = D[θ1[t], t];
ω2 = D[θ1[t] + θ2[t], t];
ω3 = D[θ1[t] + θ2[t] + θ3[t], t];
ω4 = D[θ4[t], t];
ω5 = D[Pi/2 - θ4[t] + θ5[t], t];
ω6 = D[Pi/2 - θ4[t] + θ5[t] + θ6[t], t];

T = 1/2*m1*((D[XG1, t])^2 + (D[YG1, t])^2) + 1/2*I1*ω1^2 +
1/2*m2*((D[XG2, t])^2 + (D[YG2, t])^2) + 1/2*I2*ω2^2 +
1/2*m3*((D[XG3, t])^2 + (D[YG3, t])^2) + 1/2*I3*ω3^2 +
1/2*m4*((D[XG4, t])^2 + (D[YG4, t])^2) + 1/2*I4*ω4^2 +
1/2*m5*((D[XG5, t])^2 + (D[YG5, t])^2) + 1/2*I5*ω5^2 +
1/2*m6*((D[XG6, t])^2 + (D[YG6, t])^2) + 1/2*I6*ω6^2;
V = m1*g*YG1 + m2*g*YG2 + m3*g*YG3 + m4*g*YG4 + m5*g*YG5 + m6*g*YG6;

(*Lagrangian*)
L = T - V;

(*Equations*)
(*Energy equations*)
eq1 = D[D[L, θ4'[t]], t] - D[L, θ4[t]] +
D[RDF, θ4'[t]] == T4;
eq2 = D[D[L, θ5'[t]], t] - D[L, θ5[t]] +
D[RDF, θ5'[t]] == T5;
eq3 = D[D[L, θ6'[t]], t] - D[L, θ6[t]] +
D[RDF, θ6'[t]] == T6;

ics = {θ4 == Pi/2, θ4' == 0, θ5 ==
0, θ5' == 0, θ6 == 0, θ6' == 0};
eqs = Join[{eq1, eq2, eq3}, ics];

s = NDSolve[eqs, {θ4, θ5, θ6}, {t, 0, 10},
MaxStepSize -> 0.05, AccuracyGoal -> 10, PrecisionGoal -> 10,
MaxSteps -> 1000000,
Method -> {"EquationSimplification" -> "Residual"}];

Plot[Evaluate[{θ4[t]} /. s], {t, 0, 10},
PlotRange -> {{0, 15}, {-Pi, Pi}}]
Plot[Evaluate[{θ4'[t]} /. s], {t, 0, 10},
PlotRange -> {{0, 15}, {-5, 5}}]

Plot[Evaluate[{θ5[t]} /. s], {t, 0, 10},
PlotRange -> {{0, 15}, {-5, 5}}]
Plot[Evaluate[{θ5'[t]} /. s], {t, 0, 10},
PlotRange -> {{0, 15}, {-Pi, Pi}}]

Plot[Evaluate[{θ6[t]} /. s], {t, 0, 10},
PlotRange -> {{0, 15}, {-2 Pi, 2 Pi}}]
Plot[Evaluate[{θ6'[t]} /. s], {t, 0, 10},
PlotRange -> {{0, 15}, {-Pi, Pi}}]

My dynamics should be completely stable. My equations set is unstable(I make mistake in drive equation)? Or my nonlinear solving method is not proper method?(for example i should choose a better explicit or implicit method for solve these nonlinear equations)

Also you can see that behavior(see my first plot of Theta1,Theta2 after run my codes) of system is stable but in a certain time it diverges strangly!!! • There is a typo mistake in YG3. Complete the Sin argument. – zhk Feb 12 '17 at 16:29
• I complete Sin in my code but it generate error like this: – Farshid Feb 13 '17 at 6:04
• NDSolve::mconly: For the method IDA, only machine real code is available. Unable to continue with complex values or beyond floating-point exceptions. whats the meaning – Farshid Feb 13 '17 at 6:05
• I'm voting to close this question as off-topic because it's too localized and unlikely to help future vistors. – xzczd Feb 13 '17 at 12:25
• Possibly a Sqrt of a negative number happened -- perhaps you can check. I suppose that mathematically that should not happen, in which case it might occur as a result of rounding, truncation or coding errors. But first I would want to figure out which of the two problems (complex or FP exception) is occurring. – Goofy Feb 19 '17 at 15:31