rho = 0.4;      % ratio of mass of lower piece to mass of upper piece
alpha1 = 0.4;   % normalized distances from pivot to center of mass
alpha2 = 0.4;
beta1 = 0.3;    % normalized moments of inertia
beta2 = 0.3;
g = 9.81;       % gravitational acceleration
len = 0.248;    % distance in meters between pivots

fun = @(t,y) pendrhs(t,y,[rho+alpha1 rho*alpha2 rho+beta1 rho*beta2 g/len]);
yinit = [1 0 0 0];

fps = 30;       % frames per second
secs = 100;     % number of seconds to simulate

for n = 1:fps*secs;
    [t y] = ode45(fun, [0 1/fps], yinit);
    yinit = y(end,:);
    plot([0 sin(yinit(1)) sin(yinit(1))+0.9*sin(yinit(2))],...
        [0 -cos(yinit(1)) -cos(yinit(1))-0.9*cos(yinit(2))]);
    axis([-2 2 -2 2]);
    axis square;
    pause(1/fps);
end
