(*position analysis*)

m1=m2=m;

L1=L2=L;

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]};

JC1=m1*L1^2/12.;

JC2=m2*L2^2/12.;
T1=m1 vC1.vC1/2.+JC1 w1.w1/2.;
T2=m2 vC2.vC2/2.+JC2 w2.w2/2.;

(*T1 can be calculated as T1=J0 w1.w1/2.
where  J0=m1*L1^2/3 *)
Simplify[T1];
Simplify[T2];

T=T1+T2;

Simplify[T];
(*external forces*)

FC1={0,m1 g,0};

FC2={0,m2 g,0};
(*generalized forces*)

Q1=FC1.D[rC1,q1[t]]+FC2.D[rC2,q1[t]];

Q2=FC1.D[rC1,q2[t]]+FC2.D[rC2,q2[t]];
(*Lagrange's equations*)

eq1=D[D[T,q1'[t],t]]-D[T,q1[t]]-Q1;

eq2=D[D[T,q2'[t],t]]-D[T,q2[t]]-Q2;

(*input data*)

inp={m1->1.,L1->1.,m2->1.,L2->1.,g->10.};
equation1=eq1/.inp;

equation2=eq2/.inp;
sol=NDSolve[{equation1==0,equation2==0,
q1[0]==.1,q2[0]==.1,q1'[0]==0.,q2'[0]==0.},
{q1,q2},{t,0.,40.},MaxSteps->2000];

Plot[Evaluate[q1[t]]/.sol,{t,0.,40.}]
Plot[Evaluate[q2[t]]/.sol,{t,0.,40.}]