function pend
close all
   g = 9.81;  L = 1; d = 0.1; beta = 1.85;
   Omega = beta*sqrt(g/L);
   initial_w = [0;0];
   [times,sols] = ode45(@(t,w) diffeq(t,w,Omega,g,L,d),[0:0.1:60],initial_w);
   animate_pendulum_rotor(d,Omega,L,times,sols)

end
function dwdt = diffeq(t,w,Omega,g,L,d)
   theta = w(1); omega = w(2);
   dthetadt = omega;
   domegadt = (d*Omega^2*cos(theta-Omega*t)-g*sin(theta))/L;
   dwdt = [dthetadt;domegadt];
end