function dydt = ode2(t,y)
% assume 'y' is a vector of form y = [ theta1 ; theta1_prime ; theta2 ; theta2_prime]

theta1= [y(1)];
%theta1_prime= change of theta with respect to time, omega. 
% user set omega values for simplicity.
omega1 = [y(2)];
theta2 = [y(3)];
%theta2_prime=
omega2=[y(4)];

% use the following three variables from the main worksheet
% must declare them 'global' here and there
global L1 L2 m1 m2 g


% omega_primes are derived for the equations of motion.
omega1_prime=(-m2*L1*y(2)^2*sin(-2*y(1)+2*y(3))-g*m2*sin(2*y(3)-y(1))+...
    g*sin(y(1))*m2+2*m2*L2*y(4)^2*sin(-y(1)+y(3))+2*g*sin(y(1))*m1)/...
    (-L1*m2+L1*m2*cos(-2*y(1)+2*y(3))-2*L1*m1);
omega2_prime=(g*m2*sin(y(3))+2*m2*L1*y(4)^2*sin(-y(1)+y(3))+...
    g*sin(y(3))*m1-m2*L2*y(4)^2*sin(-2*y(1)+2*y(3))+...
    2*L1*y(4)^2*m1*sin(-y(1)+y(3))+g*m1*sin(y(3)-2*y(1))+g*m2*sin(y(3)-2*y(1)))/...
    (-m2*L2+L2*m2*cos(-2*y(1)+2*y(3))-2*L2*m1);

dydt = [ omega1 ; omega1_prime ; omega2 ; omega2_prime];  % must column vector 