(* Input data *)
(* L1 = AB; L2 = BD; L3 are the lengths*)
(* h1, h2 and h3 are the widths *)
(* a is the thickness *)
(* ro[Kg/m^3] is the density of steel*)
(* m1, m2 and m3 are the masses *)
(* IC1, IC2 and IC3 are the moments of inertia *)
(* g is the gravitational acceleration *)
(* M1 is the motor torque *)

L1 = .100;
L2 = .470;
L3 = .047;
AC = .280;

h1 = h2 = h = .01;
h3 = .025;
ro = 7850;
a = .005;

m1 = ro*L1*h1*a;
m2 = ro*L2*h2*a;
m3a = ro*L3*h3*a;
m3b = ro*L3*h*a;
m3 = m3a - m3b;

IC1 = m1/12*(L1^2 + h1^2);
IC2 = m2/12*(L2^2 + h2^2);
IC3 = m3a/12*(L3^2 + h3^2) - m3b/12*(L3^2 + h^2);

g = 9.81;

(* Position, velocity and acceleration vectors *)

rC1 = {L1*Cos[theta[t]]/2., L1*Sin[theta[t]]/2., 0};
vC1 = D[rC1, t];
theta2 = ArcTan[L1*Sin[theta[t]]/(L1*Cos[theta[t]] - AC)];
rC2 = {L1*Cos[theta[t]] + L2*Cos[theta2]/2.,
      L1*Sin[theta[t]] + L2*Sin[theta2]/2., 0};
vC2 = D[rC2, t];
rC3 = {AC, 0, 0};
vC3 = D[rC3, t];

(* Angular velocities *)

omega = {0, 0, theta'[t]};
omega2 = Simplify[{0, 0, D[theta2, t]}];

(* Kinetic energy *)

T1 = m1*vC1.vC1/2. + IC1*omega.omega/2.;
T2 = m2*vC2.vC2/2. + IC2*omega2.omega2/2.;
T3 = m3*vC3.vC3/2. + IC3*omega2.omega2/2.;
T = Expand[T1 + T2 + T3];

(* Left hand side of the Lagrange's equation *)

LHS = D[D[T, theta'[t]], t] - D[T, theta[t]];

(* Right hand side of the Lagrange's equation *)(* generalized forces *)

G1 = {0, -m1  g, 0};
G2 = {0, -m2  g, 0};
G3 = {0, -m3  g, 0};

(* DC motor torque *)
M1 = {0, 0, 1. - .25  theta'[t]};

RHS = D[rC1, theta[t]].G1 + D[rC2, theta[t]].G2 + D[rC3, theta[t]].G3 +
      D[omega, theta'[t]].M1;
 

eqLHS = (LHS /. {theta'[t] -> w[t], theta''[t] -> w'[t]});
eqRHS = (RHS /. {theta'[t] -> w[t], theta''[t] -> w'[t]});

soluti = NDSolve[{eqLHS == eqRHS, theta'[t] == w[t], theta[0] == N[Pi/6],
        w[0] == 0}, {theta[t], w[t]}, {t, 0, 5}];

Plot[Evaluate[theta[t]] /. soluti, {t, 0., 5}, PlotRange -> {All, All}, AxesLabel -> {"t[s]", "theta[rad]"}];
Plot[(Evaluate[w[t]] /. soluti)* 30/N[Pi]  , {t, 0., 5}, PlotRange -> {All, {0, 55}}, AxesLabel -> {"t[s]", "n[rpm]"}];
Plot[(Evaluate[w[t]] /. soluti)* 30/N[Pi]  , {t, 0., 5}, PlotRange -> {All, {30, 43}}, AxesLabel -> {"t[s]", "n[rpm]"}];
Plot[(Evaluate[w[t]] /. soluti)* 30/N[Pi]  , {t, 0., .07},PlotRange -> {All, {0, 45}}, AxesLabel -> {"t[s]", "n[rpm]"}];