Off[General::spell1];
Off[General::spell];

(*position analysis*)

rA={L1*Sin[q1[t]],L1*Cos[q1[t]],0};

rC1={L1*Sin[q1[t]]/2.,L1*Cos[q1[t]]/2.,0};

rB=rA+{L2*Sin[q2[t]],L2*Cos[q2[t]],0};

rC2=rA+{L2*Sin[q2[t]]/2.,L2*Cos[q2[t]]/2.,0};

(*velocity analysis*)

vA=D[rA,t];

vC1=D[rC1,t];

vB=D[rB,t];

vC2=D[rC2,t];

w1={0,0,q1'[t]};

w2={0,0,q2'[t]};


(*acceleration analysis*)

aA=D[vA,t];

aC1=D[vC1,t];

aB=D[vB,t];

aC2=D[vC2,t];

alpha1={0,0,q1''[t]};

alpha2={0,0,q2''[t]};


FO={FOx,FOy,0};

FA={FAx,FAy,0};

G1={0,m1 g,0};

G2={0,m2 g,0};


JC1=m1*L1^2/12.;

JC2=m2*L2^2/12.;

N1=m1 aC1 - FO - FA - G1;
N2=m2 aC2 + FA - G2;

E1=JC1 alpha1 - Cross[(-rC1),FO] - Cross[(rA-rC1),FA];

E2=JC2 alpha2 - Cross[(rA-rC2),FA];

Simplify[N1[[1]]]
Simplify[N1[[2]]]
Simplify[N2[[1]]]
Simplify[N2[[2]]]


solu=Solve[{N1[[1]]==0,N1[[2]]==0,
  N2[[1]]==0,N2[[2]]==0},{FOx,FOy,FAx,FAy}];


(*input data*)
inp={m1->0.4,L1->1.,m2->0.8,L2->2.,g->9.8};


soluti = NDSolve[{((E1[[3]] /. solu)[[1]] /. inp) == 
        0, ((E2[[3]] /. solu)[[1]] /. inp) == 0, q1[0] == N[Pi/3], 
      q2[0] == N[Pi]/6, q1'[0] == 0., q2'[0] == 0.}, {q1, q2}, {t, 0., 10.}, 
    MaxSteps -> 2000]

solu=Solve[{N1[[1]]==0,N1[[2]]==0,
  N2[[1]]==0,N2[[2]]==0},{FOx,FOy,FAx,FAy}];

(*input data*)
inp={m1->1.,L1->1.,m2->1.,L2->1.,g->10.};


soluti = NDSolve[{
      ((E1[[3]] /. solu)[[1]] /. inp) == 0,
      ((E2[[3]] /. solu)[[1]] /. inp) == 0,
      q1[0] == N[Pi/3], q2[0] == N[Pi]/6, q1'[0] == 0., q2'[0] == 0.}, {q1, 
      q2}, {t, 0., 10.}, MaxSteps -> 2000]


Plot[Evaluate[q1[t]] /. soluti, {t, 0., 10.}]
Plot[Evaluate[q2[t]] /. soluti, {t, 0., 10.}]