function output=f_4b(theta1,theta2,omega1,omega2,t)
m1=1;
m2=2;
L1=1;
L2=2;
g=9.8;

D=L2*(2*m1+m2-m2*cos(2*theta1-2*theta2));
C=2*sin(theta1-theta2)*(omega1^2*L1*(m1+m2)+g*(m1+m2)*cos(theta1)+omega2^2*L2*m2*cos(theta1-theta2));
output=C/D;