function stated = DoublePendulum(t, state)

a1 = state(1);
a2 = state(2);
a1d = state(3);
a2d = state(4);

damping = 0.0;
g = 9.8;

% inertia matrix
M = [3 + 2*cos(a2), 1+cos(a2); ...
     1+cos(a2), 1];

% coriolis, centripetal and gravitational forces
c1 = a2d*((2*a1d) + a2d)*sin(a2) + ...
     2*g*sin(a1) + g*sin(a1+a2);
c2 = -(a1d^2)*sin(a2) + g*sin(a1+a2);

% passive dynamics
cc = [c1-damping*a1d; ...
      c2-damping*a2d];

% compute accelerations
acc = M\cc;

stated = [a1d; a2d; acc(1); acc(2)];

end
