Quadruple pendulum yo

function quadruplependulumsim
close all
t=0:0.001:60;
g=9.81;
m1=1;
m2=1;
m3=1;
m4=1;
l1=1.6;
l2=1.4;
l3=1.2;
l4=1;
theta1_initial = 0.2;
theta1dot_initial = 0;
theta2_initial = 0.4;
theta2dot_initial = 0;
theta3_initial = 0.6;
theta3dot_initial = 0;
theta4_initial = 0.8;
theta4dot_initial = 0;
[t,x]=ode45(@rhs,t,[theta1_initial theta1dot_initial theta2_initial theta2dot_initial theta3_initial theta3dot_initial theta4_initial theta4dot_initial]);
plot(t,x(:,1),t,x(:,3),t,x(:,5),t,x(:,7));
xlabel('t (s)'),ylabel('Angle (rad)');
%pause
    function f = rhs(t,x)
        f_1 = x(2);
        f_2 = -((g*m2^2*m3*sin(2*x(3) - x(1) + 2*x(5)))/4 - (g*m2^2*m3*sin(x(1) - 2*x(3)))/2 - (g*m2^2*m4*sin(x(1) - 2*x(3)))/4 - (g*m2*m3^2*sin(x(1) - 2*x(5)))/2 - (g*m2^2*m3*sin(x(1) - 2*x(5)))/4 - (g*m2^2*m3*sin(x(1) + 2*x(5)))/4 - (g*m2*m3^2*sin(x(1) - 2*x(3)))/2 + (g*m2^2*m4*sin(x(1) - 2*x(3) - 2*x(5) + 2*x(7)))/8 + (g*m2^2*m4*sin(x(1) - 2*x(3) + 2*x(5) - 2*x(7)))/8 - (g*m2*m3^2*sin(x(1)))/2 - (g*m2^2*m3*sin(x(1)))/2 - (g*m2^2*m4*sin(x(1)))/4 - (g*m2*m3^2*sin(x(1) - 2*x(3) + 2*x(5)))/2 - (g*m2^2*m3*sin(x(1) - 2*x(3) + 2*x(5)))/4 + (g*m2^2*m4*sin(x(1) - 2*x(5) + 2*x(7)))/8 + (g*m2^2*m4*sin(x(1) + 2*x(5) - 2*x(7)))/8 + (g*m1*m3*m4*sin(x(1) - 2*x(3) + 2*x(5)))/4 + (g*m1*m3*m4*sin(x(1) + 2*x(3) - 2*x(5)))/4 - (g*m2*m3*m4*sin(x(1) - 2*x(3) + 2*x(5)))/8 + (g*m1*m2*m4*sin(x(1) - 2*x(5) + 2*x(7)))/4 + (g*m1*m2*m4*sin(x(1) + 2*x(5) - 2*x(7)))/4 - (g*m2*m3*m4*sin(x(1) - 2*x(3) + 2*x(7)))/8 + (g*m2*m3*m4*sin(x(1) - 2*x(5) + 2*x(7)))/8 - (l1*m2^2*m4*sin(2*x(1) - 2*x(3) - 2*x(5) + 2*x(7))*x(2)^2)/8 - (l1*m2^2*m4*sin(2*x(1) - 2*x(3) + 2*x(5) - 2*x(7))*x(2)^2)/8 - (g*m1*m2*m3*sin(x(1) - 2*x(5)))/2 - (g*m1*m2*m3*sin(x(1) + 2*x(5)))/2 - (5*g*m2*m3*m4*sin(x(1) - 2*x(3)))/8 - (g*m2*m3*m4*sin(x(1) - 2*x(5)))/8 - (g*m1*m3*m4*sin(x(1) - 2*x(7)))/4 - (g*m1*m3*m4*sin(x(1) + 2*x(7)))/4 - (g*m2*m3*m4*sin(x(1) - 2*x(7)))/8 - (g*m2*m3*m4*sin(x(1) + 2*x(7)))/4 - (l1*m2^2*m3*sin(2*x(3) - 2*x(1) + 2*x(5))*x(2)^2)/4 + (l1*m2^2*m3*sin(2*x(1) - 2*x(3) + 2*x(5))*x(2)^2)/4 + (l1*m2^2*m3*x(2)^2*sin(2*x(1) - 2*x(3)))/2 + (l1*m2^2*m4*x(2)^2*sin(2*x(1) - 2*x(3)))/4 + (g*m2*m3*m4*sin(2*x(3) - x(1) + 2*x(7)))/4 - (l2*m2^2*m4*x(4)^2*sin(x(1) - x(3) - 2*x(5) + 2*x(7)))/4 - (l2*m2^2*m4*x(4)^2*sin(x(1) - x(3) + 2*x(5) - 2*x(7)))/4 + (l2*m2^2*m3*x(4)^2*sin(x(1) - x(3) + 2*x(5)))/2 - (l2*m2^2*m3*x(4)^2*sin(x(3) - x(1) + 2*x(5)))/2 + (g*m1*m3*m4*sin(x(1) - 2*x(3) + 2*x(5) - 2*x(7)))/8 + (g*m1*m3*m4*sin(x(1) - 2*x(3) + 2*x(5) + 2*x(7)))/8 + (g*m1*m3*m4*sin(x(1) + 2*x(3) - 2*x(5) - 2*x(7)))/8 + (g*m1*m3*m4*sin(x(1) + 2*x(3) - 2*x(5) + 2*x(7)))/8 + (g*m2*m3*m4*sin(x(1) - 2*x(3) + 2*x(5) - 2*x(7)))/8 - g*m1*m2*m3*sin(x(1)) - (g*m1*m2*m4*sin(x(1)))/2 - (g*m1*m3*m4*sin(x(1)))/2 - (5*g*m2*m3*m4*sin(x(1)))/8 + l2*m2^2*m3*sin(x(1) - x(3))*x(4)^2 + (l2*m2^2*m4*sin(x(1) - x(3))*x(4)^2)/2 + l2*m2*m3*m4*sin(x(1) - x(3))*x(4)^2 + (l3*m2*m3*m4*sin(x(1) - x(5))*x(6)^2)/2 + (l4*m2*m3*m4*sin(x(1) - x(7))*x(8)^2)/4 - (l1*m2*m3*m4*sin(2*x(3) - 2*x(1) + 2*x(7))*x(2)^2)/4 + (l1*m2*m3*m4*sin(2*x(1) - 2*x(3) + 2*x(7))*x(2)^2)/4 + (l3*m2*m3*m4*sin(x(1) - 2*x(3) + x(5) - 2*x(7))*x(6)^2)/4 + (l3*m2*m3*m4*sin(x(1) - 2*x(3) + x(5) + 2*x(7))*x(6)^2)/4 + (l4*m2*m3*m4*sin(x(1) - 2*x(3) + 2*x(5) + x(7))*x(8)^2)/4 + (l1*m2*m3*m4*x(2)^2*sin(2*x(1) - 2*x(3)))/2 + (l3*m2*m3*m4*x(6)^2*sin(x(1) - 2*x(3) + x(5)))/2 + (l4*m2*m3*m4*x(8)^2*sin(x(1) - 2*x(3) + x(7)))/4 + (l4*m2*m3*m4*x(8)^2*sin(x(1) - 2*x(5) + x(7)))/4 + (l4*m2*m3*m4*x(8)^2*sin(x(1) + x(7)))/4 + (l4*m2*m3*m4*x(8)^2*sin(x(1) - 2*x(3) + 2*x(5) - x(7)))/4 + (l2*m2*m3*m4*x(4)^2*sin(x(1) - x(3) + 2*x(7)))/2 - (l2*m2*m3*m4*x(4)^2*sin(x(3) - x(1) + 2*x(7)))/2 + (l3*m2*m3*m4*x(6)^2*sin(x(1) - x(5) + 2*x(7)))/4 - (l3*m2*m3*m4*x(6)^2*sin(x(5) - x(1) + 2*x(7)))/4 - (l4*m2*m3*m4*x(8)^2*sin(2*x(3) - x(1) + x(7)))/4 - (l4*m2*m3*m4*x(8)^2*sin(2*x(5) - x(1) + x(7)))/4)/(l1*((m2^2*m3)/2 + (m2^2*m4)/4 - (m2^2*m3*cos(2*x(1) - 2*x(3)))/2 - (m2^2*m4*cos(2*x(1) - 2*x(3)))/4 - (m2^2*m4*cos(2*x(5) - 2*x(7)))/4 + m1*m2*m3 + (m1*m2*m4)/2 + (m1*m3*m4)/2 + (m2*m3*m4)/2 + (m2^2*m4*cos(2*x(1) - 2*x(3) - 2*x(5) + 2*x(7)))/8 + (m2^2*m4*cos(2*x(1) - 2*x(3) + 2*x(5) - 2*x(7)))/8 + (m2^2*m3*cos(2*x(5)))/2 - (m2^2*m3*cos(2*x(3) - 2*x(1) + 2*x(5)))/4 - (m2^2*m3*cos(2*x(1) - 2*x(3) + 2*x(5)))/4 - (m2*m3*m4*cos(2*x(1) - 2*x(3)))/2 - (m1*m3*m4*cos(2*x(3) - 2*x(5)))/2 - (m1*m2*m4*cos(2*x(5) - 2*x(7)))/2 + m1*m2*m3*cos(2*x(5)) + (m1*m3*m4*cos(2*x(7)))/2 + (m2*m3*m4*cos(2*x(7)))/2 - (m2*m3*m4*cos(2*x(3) - 2*x(1) + 2*x(7)))/4 - (m2*m3*m4*cos(2*x(1) - 2*x(3) + 2*x(7)))/4 - (m1*m3*m4*cos(2*x(5) - 2*x(3) + 2*x(7)))/4 - (m1*m3*m4*cos(2*x(3) - 2*x(5) + 2*x(7)))/4));
        f_3 = x(4);
        f_4 = (g*m1*m3^2*sin(x(3) - 2*x(5)) + (g*m2*m3^2*sin(x(3) - 2*x(5)))/2 + (g*m2^2*m3*sin(x(3) - 2*x(5)))/4 + (g*m2^2*m3*sin(x(3) + 2*x(5)))/4 + (g*m2^2*m4*sin(2*x(1) - x(3) - 2*x(5) + 2*x(7)))/8 + (g*m2^2*m4*sin(2*x(1) - x(3) + 2*x(5) - 2*x(7)))/8 - (g*m2^2*m3*sin(2*x(1) - x(3) + 2*x(5)))/4 - (g*m2*m3^2*sin(2*x(1) - x(3)))/2 - (g*m2^2*m3*sin(2*x(1) - x(3)))/2 - (g*m2^2*m4*sin(2*x(1) - x(3)))/4 + g*m1*m3^2*sin(x(3)) + (g*m2*m3^2*sin(x(3)))/2 + (g*m2^2*m3*sin(x(3)))/2 + (g*m2^2*m4*sin(x(3)))/4 + (g*m2*m3^2*sin(x(3) - 2*x(1) + 2*x(5)))/2 + (g*m2^2*m3*sin(x(3) - 2*x(1) + 2*x(5)))/4 - (g*m2^2*m4*sin(x(3) - 2*x(5) + 2*x(7)))/8 - (g*m2^2*m4*sin(x(3) + 2*x(5) - 2*x(7)))/8 + (g*m1*m2*m3*sin(x(3) - 2*x(1) + 2*x(5)))/4 + (g*m1*m3*m4*sin(2*x(1) + x(3) - 2*x(5)))/4 + (g*m2*m3*m4*sin(x(3) - 2*x(1) + 2*x(5)))/8 + (g*m1*m3*m4*sin(x(3) - 2*x(1) + 2*x(7)))/8 + (g*m2*m3*m4*sin(x(3) - 2*x(1) + 2*x(7)))/8 - (g*m1*m2*m4*sin(x(3) - 2*x(5) + 2*x(7)))/8 - (g*m1*m2*m4*sin(x(3) + 2*x(5) - 2*x(7)))/8 - (g*m1*m3*m4*sin(x(3) - 2*x(5) + 2*x(7)))/8 - (g*m2*m3*m4*sin(x(3) - 2*x(5) + 2*x(7)))/8 - (l2*m2^2*m4*sin(2*x(1) - 2*x(3) - 2*x(5) + 2*x(7))*x(4)^2)/8 - (l2*m2^2*m4*sin(2*x(1) - 2*x(3) + 2*x(5) - 2*x(7))*x(4)^2)/8 + (g*m1*m2*m3*sin(x(3) - 2*x(5)))/4 + (g*m1*m2*m3*sin(x(3) + 2*x(5)))/4 + (g*m1*m3*m4*sin(x(3) - 2*x(5)))/2 + (g*m2*m3*m4*sin(x(3) - 2*x(5)))/8 - (g*m1*m3*m4*sin(x(3) - 2*x(7)))/8 + (g*m1*m3*m4*sin(x(3) + 2*x(7)))/8 + (g*m2*m3*m4*sin(x(3) - 2*x(7)))/8 + (g*m2*m3*m4*sin(x(3) + 2*x(7)))/4 - (l2*m2^2*m3*sin(2*x(3) - 2*x(1) + 2*x(5))*x(4)^2)/4 + (l2*m2^2*m3*sin(2*x(1) - 2*x(3) + 2*x(5))*x(4)^2)/4 + (g*m1*m2*m4*sin(2*x(1) - x(3) - 2*x(5) + 2*x(7)))/8 + (g*m1*m2*m4*sin(2*x(1) - x(3) + 2*x(5) - 2*x(7)))/8 + (g*m2*m3*m4*sin(2*x(1) - x(3) - 2*x(5) + 2*x(7)))/8 + (l2*m2^2*m3*x(4)^2*sin(2*x(1) - 2*x(3)))/2 + (l2*m2^2*m4*x(4)^2*sin(2*x(1) - 2*x(3)))/4 - (g*m1*m2*m3*sin(2*x(1) - x(3) + 2*x(5)))/4 - (g*m1*m3*m4*sin(2*x(1) - x(3) + 2*x(7)))/8 - (g*m2*m3*m4*sin(2*x(1) - x(3) + 2*x(7)))/4 - (g*m1*m3*m4*sin(2*x(5) - x(3) + 2*x(7)))/8 - (g*m1*m2*m3*sin(2*x(1) - x(3)))/2 - (g*m1*m2*m4*sin(2*x(1) - x(3)))/4 - (g*m1*m3*m4*sin(2*x(1) - x(3)))/4 - (5*g*m2*m3*m4*sin(2*x(1) - x(3)))/8 - (l1*m2^2*m4*x(2)^2*sin(x(1) - x(3) - 2*x(5) + 2*x(7)))/4 - (l1*m2^2*m4*x(2)^2*sin(x(1) - x(3) + 2*x(5) - 2*x(7)))/4 + (l1*m2^2*m3*x(2)^2*sin(x(1) - x(3) + 2*x(5)))/2 - (l1*m2^2*m3*x(2)^2*sin(x(3) - x(1) + 2*x(5)))/2 + (g*m1*m3*m4*sin(2*x(1) + x(3) - 2*x(5) - 2*x(7)))/8 + (g*m1*m3*m4*sin(2*x(1) + x(3) - 2*x(5) + 2*x(7)))/8 + (g*m1*m2*m3*sin(x(3)))/2 + (g*m1*m2*m4*sin(x(3)))/4 + (g*m1*m3*m4*sin(x(3)))/2 + (5*g*m2*m3*m4*sin(x(3)))/8 + l1*m2^2*m3*sin(x(1) - x(3))*x(2)^2 + (l1*m2^2*m4*sin(x(1) - x(3))*x(2)^2)/2 + l1*m1*m2*m3*sin(x(1) - x(3))*x(2)^2 + (l1*m1*m2*m4*sin(x(1) - x(3))*x(2)^2)/2 + (l1*m1*m3*m4*sin(x(1) - x(3))*x(2)^2)/2 + l1*m2*m3*m4*sin(x(1) - x(3))*x(2)^2 - l3*m1*m3*m4*sin(x(3) - x(5))*x(6)^2 - (l3*m2*m3*m4*sin(x(3) - x(5))*x(6)^2)/2 - (l4*m1*m3*m4*sin(x(3) - x(7))*x(8)^2)/2 - (l4*m2*m3*m4*sin(x(3) - x(7))*x(8)^2)/4 + (l3*m2*m3*m4*sin(2*x(1) - x(3) - x(5) + 2*x(7))*x(6)^2)/4 - (l2*m2*m3*m4*sin(2*x(3) - 2*x(1) + 2*x(7))*x(4)^2)/4 + (l2*m2*m3*m4*sin(2*x(1) - 2*x(3) + 2*x(7))*x(4)^2)/4 + (l2*m1*m3*m4*sin(2*x(5) - 2*x(3) + 2*x(7))*x(4)^2)/4 - (l2*m1*m3*m4*sin(2*x(3) - 2*x(5) + 2*x(7))*x(4)^2)/4 - (l1*m1*m3*m4*sin(x(1) + x(3) - 2*x(5) - 2*x(7))*x(2)^2)/4 - (l1*m1*m3*m4*sin(x(1) + x(3) - 2*x(5) + 2*x(7))*x(2)^2)/4 - (l3*m2*m3*m4*sin(x(3) - 2*x(1) + x(5) + 2*x(7))*x(6)^2)/4 - (l4*m2*m3*m4*sin(x(3) - 2*x(1) + 2*x(5) + x(7))*x(8)^2)/4 + (l2*m2*m3*m4*x(4)^2*sin(2*x(1) - 2*x(3)))/2 - (l2*m1*m3*m4*x(4)^2*sin(2*x(3) - 2*x(5)))/2 - (l1*m1*m3*m4*x(2)^2*sin(x(1) + x(3) - 2*x(5)))/2 - (l3*m2*m3*m4*x(6)^2*sin(x(3) - 2*x(1) + x(5)))/2 - (l4*m2*m3*m4*x(8)^2*sin(x(3) - 2*x(1) + x(7)))/4 - (l4*m1*m3*m4*x(8)^2*sin(x(3) - 2*x(5) + x(7)))/2 - (l4*m2*m3*m4*x(8)^2*sin(x(3) - 2*x(5) + x(7)))/4 - (l4*m1*m3*m4*x(8)^2*sin(x(3) + x(7)))/2 - (l4*m2*m3*m4*x(8)^2*sin(x(3) + x(7)))/4 - (l1*m1*m2*m4*x(2)^2*sin(x(1) - x(3) - 2*x(5) + 2*x(7)))/4 - (l1*m1*m2*m4*x(2)^2*sin(x(1) - x(3) + 2*x(5) - 2*x(7)))/4 + (l4*m2*m3*m4*x(8)^2*sin(2*x(1) - x(3) - 2*x(5) + x(7)))/4 + (l1*m1*m2*m3*x(2)^2*sin(x(1) - x(3) + 2*x(5)))/2 - (l1*m1*m2*m3*x(2)^2*sin(x(3) - x(1) + 2*x(5)))/2 + (l1*m1*m3*m4*x(2)^2*sin(x(1) - x(3) + 2*x(7)))/4 - (l1*m1*m3*m4*x(2)^2*sin(x(3) - x(1) + 2*x(7)))/4 + (l1*m2*m3*m4*x(2)^2*sin(x(1) - x(3) + 2*x(7)))/2 - (l1*m2*m3*m4*x(2)^2*sin(x(3) - x(1) + 2*x(7)))/2 - (l3*m1*m3*m4*x(6)^2*sin(x(3) - x(5) + 2*x(7)))/2 + (l3*m1*m3*m4*x(6)^2*sin(x(5) - x(3) + 2*x(7)))/2 + (l4*m2*m3*m4*x(8)^2*sin(2*x(1) - x(3) + x(7)))/4 - (l3*m2*m3*m4*x(6)^2*sin(x(3) - x(5) + 2*x(7)))/4 + (l3*m2*m3*m4*x(6)^2*sin(x(5) - x(3) + 2*x(7)))/4 + (l4*m1*m3*m4*x(8)^2*sin(2*x(5) - x(3) + x(7)))/2 + (l4*m2*m3*m4*x(8)^2*sin(2*x(5) - x(3) + x(7)))/4)/(l2*((m2^2*m3)/2 + (m2^2*m4)/4 - (m2^2*m3*cos(2*x(1) - 2*x(3)))/2 - (m2^2*m4*cos(2*x(1) - 2*x(3)))/4 - (m2^2*m4*cos(2*x(5) - 2*x(7)))/4 + m1*m2*m3 + (m1*m2*m4)/2 + (m1*m3*m4)/2 + (m2*m3*m4)/2 + (m2^2*m4*cos(2*x(1) - 2*x(3) - 2*x(5) + 2*x(7)))/8 + (m2^2*m4*cos(2*x(1) - 2*x(3) + 2*x(5) - 2*x(7)))/8 + (m2^2*m3*cos(2*x(5)))/2 - (m2^2*m3*cos(2*x(3) - 2*x(1) + 2*x(5)))/4 - (m2^2*m3*cos(2*x(1) - 2*x(3) + 2*x(5)))/4 - (m2*m3*m4*cos(2*x(1) - 2*x(3)))/2 - (m1*m3*m4*cos(2*x(3) - 2*x(5)))/2 - (m1*m2*m4*cos(2*x(5) - 2*x(7)))/2 + m1*m2*m3*cos(2*x(5)) + (m1*m3*m4*cos(2*x(7)))/2 + (m2*m3*m4*cos(2*x(7)))/2 - (m2*m3*m4*cos(2*x(3) - 2*x(1) + 2*x(7)))/4 - (m2*m3*m4*cos(2*x(1) - 2*x(3) + 2*x(7)))/4 - (m1*m3*m4*cos(2*x(5) - 2*x(3) + 2*x(7)))/4 - (m1*m3*m4*cos(2*x(3) - 2*x(5) + 2*x(7)))/4));
        f_5 = x(6);
        f_6 = -(g*m1*m3^2*sin(2*x(3) - x(5)) - g*m1*m3^2*sin(x(5)) - g*m2*m3^2*sin(x(5)) - (g*m2^2*m3*sin(x(5)))/2 + (g*m2*m3^2*sin(2*x(3) - 2*x(1) + x(5)))/2 + (g*m2*m3^2*sin(2*x(1) - 2*x(3) + x(5)))/2 + (g*m2^2*m3*sin(2*x(3) - 2*x(1) + x(5)))/4 + (g*m2^2*m3*sin(2*x(1) - 2*x(3) + x(5)))/4 + (g*m1*m2*m3*sin(2*x(3) - 2*x(1) + x(5)))/4 - (g*m1*m2*m3*sin(2*x(1) - 2*x(3) + x(5)))/4 - (g*m1*m2*m4*sin(2*x(1) - 2*x(3) + x(5)))/8 - (g*m1*m3*m4*sin(2*x(1) - 2*x(3) + x(5)))/4 + (g*m2*m3*m4*sin(2*x(3) - 2*x(1) + x(5)))/8 + (g*m2*m3*m4*sin(2*x(1) - 2*x(3) + x(5)))/8 - (g*m1*m2*m4*sin(2*x(1) + x(5) - 2*x(7)))/8 - (g*m1*m3*m4*sin(x(5) - 2*x(1) + 2*x(7)))/8 - (g*m1*m2*m4*sin(2*x(3) + x(5) - 2*x(7)))/8 + (g*m1*m3*m4*sin(x(5) - 2*x(3) + 2*x(7)))/8 + (l3*m2^2*m4*sin(2*x(1) - 2*x(3) - 2*x(5) + 2*x(7))*x(6)^2)/8 - (l3*m2^2*m4*sin(2*x(1) - 2*x(3) + 2*x(5) - 2*x(7))*x(6)^2)/8 - (l3*m2^2*m3*sin(2*x(5))*x(6)^2)/2 + (g*m1*m2*m3*sin(2*x(1) + x(5)))/4 + (g*m1*m2*m3*sin(2*x(3) + x(5)))/4 - (g*m1*m2*m4*sin(x(5) - 2*x(7)))/8 + (g*m1*m3*m4*sin(x(5) - 2*x(7)))/8 - (g*m1*m3*m4*sin(x(5) + 2*x(7)))/8 + (g*m2*m3*m4*sin(x(5) - 2*x(7)))/4 + (l3*m2^2*m3*sin(2*x(3) - 2*x(1) + 2*x(5))*x(6)^2)/4 + (l3*m2^2*m3*sin(2*x(1) - 2*x(3) + 2*x(5))*x(6)^2)/4 + (g*m1*m2*m4*sin(2*x(1) - 2*x(3) - x(5) + 2*x(7)))/8 + (g*m2*m3*m4*sin(2*x(1) - 2*x(3) - x(5) + 2*x(7)))/8 + (l3*m2^2*m4*x(6)^2*sin(2*x(5) - 2*x(7)))/4 + (g*m1*m3*m4*sin(2*x(1) - x(5) + 2*x(7)))/8 + (g*m1*m3*m4*sin(2*x(3) - x(5) + 2*x(7)))/8 + (g*m1*m2*m3*sin(2*x(1) - x(5)))/4 + (g*m1*m2*m4*sin(2*x(1) - x(5)))/8 + (g*m1*m2*m3*sin(2*x(3) - x(5)))/4 + (g*m1*m3*m4*sin(2*x(1) - x(5)))/4 + (g*m1*m2*m4*sin(2*x(3) - x(5)))/8 + (g*m1*m3*m4*sin(2*x(3) - x(5)))/2 - (l4*m2^2*m4*x(8)^2*sin(2*x(1) - 2*x(3) + x(5) - x(7)))/4 + (l4*m2^2*m4*x(8)^2*sin(2*x(1) - 2*x(3) - x(5) + x(7)))/4 - (g*m1*m3*m4*sin(2*x(1) - 2*x(3) + x(5) - 2*x(7)))/8 - (g*m1*m3*m4*sin(2*x(1) - 2*x(3) + x(5) + 2*x(7)))/8 - (g*m2*m3*m4*sin(2*x(1) - 2*x(3) + x(5) - 2*x(7)))/8 - g*m1*m2*m3*sin(x(5)) - (g*m1*m2*m4*sin(x(5)))/8 - (g*m1*m3*m4*sin(x(5)))/2 - (g*m2*m3*m4*sin(x(5)))/4 + (l4*m2^2*m4*sin(x(5) - x(7))*x(8)^2)/2 - (l1*m1*m2*m3*sin(x(1) - x(5))*x(2)^2)/2 - (l1*m1*m2*m4*sin(x(1) - x(5))*x(2)^2)/4 - (l1*m1*m3*m4*sin(x(1) - x(5))*x(2)^2)/2 - l2*m1*m2*m3*sin(x(3) - x(5))*x(4)^2 - (l2*m1*m2*m4*sin(x(3) - x(5))*x(4)^2)/2 - l2*m1*m3*m4*sin(x(3) - x(5))*x(4)^2 + l4*m1*m2*m4*sin(x(5) - x(7))*x(8)^2 + (l4*m1*m3*m4*sin(x(5) - x(7))*x(8)^2)/2 + (l4*m2*m3*m4*sin(x(5) - x(7))*x(8)^2)/2 - l3*m1*m2*m3*sin(2*x(5))*x(6)^2 + (l3*m1*m3*m4*sin(2*x(5) - 2*x(3) + 2*x(7))*x(6)^2)/4 - (l3*m1*m3*m4*sin(2*x(3) - 2*x(5) + 2*x(7))*x(6)^2)/4 + (l1*m1*m3*m4*sin(x(1) - 2*x(3) + x(5) - 2*x(7))*x(2)^2)/4 + (l1*m1*m3*m4*sin(x(1) - 2*x(3) + x(5) + 2*x(7))*x(2)^2)/4 - (l4*m2*m3*m4*sin(2*x(3) - 2*x(1) + x(5) + x(7))*x(8)^2)/4 - (l4*m2*m3*m4*sin(2*x(1) - 2*x(3) + x(5) + x(7))*x(8)^2)/4 - (l3*m1*m3*m4*x(6)^2*sin(2*x(3) - 2*x(5)))/2 + (l3*m1*m2*m4*x(6)^2*sin(2*x(5) - 2*x(7)))/2 + (l1*m1*m2*m3*x(2)^2*sin(x(1) - 2*x(3) + x(5)))/2 + (l1*m1*m2*m4*x(2)^2*sin(x(1) - 2*x(3) + x(5)))/4 + (l1*m1*m3*m4*x(2)^2*sin(x(1) - 2*x(3) + x(5)))/2 + (l1*m1*m2*m4*x(2)^2*sin(x(1) + x(5) - 2*x(7)))/4 + (l2*m1*m2*m4*x(4)^2*sin(x(3) + x(5) - 2*x(7)))/2 + (l4*m1*m3*m4*x(8)^2*sin(x(5) - 2*x(3) + x(7)))/2 - (l1*m1*m2*m3*x(2)^2*sin(x(1) + x(5)))/2 - l2*m1*m2*m3*x(4)^2*sin(x(3) + x(5)) + (l4*m1*m3*m4*x(8)^2*sin(x(5) + x(7)))/2 + (l4*m2*m3*m4*x(8)^2*sin(x(5) + x(7)))/2 - (l1*m1*m2*m4*x(2)^2*sin(x(1) - 2*x(3) - x(5) + 2*x(7)))/4 - (l4*m2*m3*m4*x(8)^2*sin(2*x(1) - 2*x(3) + x(5) - x(7)))/4 + (l4*m2*m3*m4*x(8)^2*sin(2*x(1) - 2*x(3) - x(5) + x(7)))/4 - (l1*m1*m2*m3*x(2)^2*sin(2*x(3) - x(1) + x(5)))/2 - (l1*m1*m3*m4*x(2)^2*sin(x(1) - x(5) + 2*x(7)))/4 + (l1*m1*m3*m4*x(2)^2*sin(x(5) - x(1) + 2*x(7)))/4 - (l2*m1*m3*m4*x(4)^2*sin(x(3) - x(5) + 2*x(7)))/2 + (l2*m1*m3*m4*x(4)^2*sin(x(5) - x(3) + 2*x(7)))/2 - (l4*m1*m3*m4*x(8)^2*sin(2*x(3) - x(5) + x(7)))/2)/(l3*((m2^2*m3)/2 + (m2^2*m4)/4 - (m2^2*m3*cos(2*x(1) - 2*x(3)))/2 - (m2^2*m4*cos(2*x(1) - 2*x(3)))/4 - (m2^2*m4*cos(2*x(5) - 2*x(7)))/4 + m1*m2*m3 + (m1*m2*m4)/2 + (m1*m3*m4)/2 + (m2*m3*m4)/2 + (m2^2*m4*cos(2*x(1) - 2*x(3) - 2*x(5) + 2*x(7)))/8 + (m2^2*m4*cos(2*x(1) - 2*x(3) + 2*x(5) - 2*x(7)))/8 + (m2^2*m3*cos(2*x(5)))/2 - (m2^2*m3*cos(2*x(3) - 2*x(1) + 2*x(5)))/4 - (m2^2*m3*cos(2*x(1) - 2*x(3) + 2*x(5)))/4 - (m2*m3*m4*cos(2*x(1) - 2*x(3)))/2 - (m1*m3*m4*cos(2*x(3) - 2*x(5)))/2 - (m1*m2*m4*cos(2*x(5) - 2*x(7)))/2 + m1*m2*m3*cos(2*x(5)) + (m1*m3*m4*cos(2*x(7)))/2 + (m2*m3*m4*cos(2*x(7)))/2 - (m2*m3*m4*cos(2*x(3) - 2*x(1) + 2*x(7)))/4 - (m2*m3*m4*cos(2*x(1) - 2*x(3) + 2*x(7)))/4 - (m1*m3*m4*cos(2*x(5) - 2*x(3) + 2*x(7)))/4 - (m1*m3*m4*cos(2*x(3) - 2*x(5) + 2*x(7)))/4));
        f_7 = x(8);
        f_8 = (g*m2^2*m3*sin(2*x(1) - 2*x(3) + 2*x(5) - x(7)) - 2*g*m2^2*m3*sin(2*x(5) - x(7)) - g*m2^2*m3*sin(2*x(1) - 2*x(3) - 2*x(5) + x(7)) - 8*g*m1*m3^2*sin(x(7)) - 8*g*m2*m3^2*sin(x(7)) - 2*g*m2^2*m3*sin(x(7)) + 4*g*m2*m3^2*sin(2*x(3) - 2*x(1) + x(7)) + 4*g*m2*m3^2*sin(2*x(1) - 2*x(3) + x(7)) + g*m2^2*m3*sin(2*x(3) - 2*x(1) + x(7)) + g*m2^2*m3*sin(2*x(1) - 2*x(3) + x(7)) + 4*g*m1*m3^2*sin(2*x(5) - 2*x(3) + x(7)) + 4*g*m1*m3^2*sin(2*x(3) - 2*x(5) + x(7)) + g*m1*m2*m3*sin(2*x(3) - 2*x(1) + x(7)) + g*m1*m2*m3*sin(2*x(1) - 2*x(3) + x(7)) + g*m1*m2*m4*sin(2*x(1) - 2*x(3) + x(7)) + g*m1*m2*m3*sin(2*x(5) - 2*x(1) + x(7)) + g*m1*m2*m3*sin(2*x(1) - 2*x(5) + x(7)) + g*m1*m2*m4*sin(2*x(1) - 2*x(5) + x(7)) + 2*g*m2*m3*m4*sin(2*x(3) - 2*x(1) + x(7)) + 2*g*m2*m3*m4*sin(2*x(1) - 2*x(3) + x(7)) + g*m1*m2*m3*sin(2*x(5) - 2*x(3) + x(7)) + g*m1*m2*m3*sin(2*x(3) - 2*x(5) + x(7)) + g*m1*m2*m4*sin(2*x(3) - 2*x(5) + x(7)) + 2*g*m1*m3*m4*sin(2*x(5) - 2*x(3) + x(7)) + 2*g*m1*m3*m4*sin(2*x(3) - 2*x(5) + x(7)) + l4*m2^2*m4*sin(2*x(1) - 2*x(3) - 2*x(5) + 2*x(7))*x(8)^2 - l4*m2^2*m4*sin(2*x(1) - 2*x(3) + 2*x(5) - 2*x(7))*x(8)^2 + g*m1*m2*m3*sin(2*x(1) + x(7)) + g*m1*m2*m3*sin(2*x(3) + x(7)) + g*m1*m2*m3*sin(2*x(5) + x(7)) - g*m1*m2*m3*sin(2*x(1) - 2*x(3) + 2*x(5) - x(7)) - g*m1*m2*m4*sin(2*x(1) - 2*x(3) + 2*x(5) - x(7)) + 2*l3*m2^2*m3*sin(2*x(3) - 2*x(1) + x(5) + x(7))*x(6)^2 + 2*l3*m2^2*m3*sin(2*x(1) - 2*x(3) + x(5) + x(7))*x(6)^2 + 2*l4*m2^2*m4*x(8)^2*sin(2*x(5) - 2*x(7)) - g*m1*m2*m3*sin(2*x(1) - x(7)) - g*m1*m2*m4*sin(2*x(1) - x(7)) - g*m1*m2*m3*sin(2*x(3) - x(7)) - g*m1*m2*m4*sin(2*x(3) - x(7)) - 5*g*m1*m2*m3*sin(2*x(5) - x(7)) - g*m1*m2*m4*sin(2*x(5) - x(7)) - 4*l3*m2^2*m3*x(6)^2*sin(x(5) + x(7)) - 2*l3*m2^2*m3*x(6)^2*sin(2*x(1) - 2*x(3) + x(5) - x(7)) + 2*l3*m2^2*m3*x(6)^2*sin(2*x(1) - 2*x(3) - x(5) + x(7)) - 2*l3*m2^2*m4*x(6)^2*sin(2*x(1) - 2*x(3) + x(5) - x(7)) + 2*l3*m2^2*m4*x(6)^2*sin(2*x(1) - 2*x(3) - x(5) + x(7)) + g*m1*m2*m3*sin(2*x(1) - 2*x(3) + 2*x(5) + x(7)) - 2*g*m1*m2*m3*sin(x(7)) + g*m1*m2*m4*sin(x(7)) - 4*g*m1*m3*m4*sin(x(7)) - 4*g*m2*m3*m4*sin(x(7)) + 4*l3*m2^2*m3*sin(x(5) - x(7))*x(6)^2 + 4*l3*m2^2*m4*sin(x(5) - x(7))*x(6)^2 + 2*l1*m1*m2*m3*sin(x(1) - x(7))*x(2)^2 + 2*l1*m1*m2*m4*sin(x(1) - x(7))*x(2)^2 + 4*l2*m1*m2*m3*sin(x(3) - x(7))*x(4)^2 + 4*l2*m1*m2*m4*sin(x(3) - x(7))*x(4)^2 + 8*l3*m1*m2*m3*sin(x(5) - x(7))*x(6)^2 + 8*l3*m1*m2*m4*sin(x(5) - x(7))*x(6)^2 + 4*l4*m1*m3*m4*sin(2*x(7))*x(8)^2 + 4*l4*m2*m3*m4*sin(2*x(7))*x(8)^2 - 2*l4*m2*m3*m4*sin(2*x(3) - 2*x(1) + 2*x(7))*x(8)^2 - 2*l4*m2*m3*m4*sin(2*x(1) - 2*x(3) + 2*x(7))*x(8)^2 - 2*l4*m1*m3*m4*sin(2*x(5) - 2*x(3) + 2*x(7))*x(8)^2 - 2*l4*m1*m3*m4*sin(2*x(3) - 2*x(5) + 2*x(7))*x(8)^2 - 2*l1*m1*m2*m3*sin(x(1) - 2*x(3) + 2*x(5) + x(7))*x(2)^2 + 4*l4*m1*m2*m4*x(8)^2*sin(2*x(5) - 2*x(7)) - 2*l1*m1*m2*m3*x(2)^2*sin(x(1) - 2*x(3) + x(7)) - 2*l1*m1*m2*m4*x(2)^2*sin(x(1) - 2*x(3) + x(7)) - 2*l1*m1*m2*m3*x(2)^2*sin(x(1) - 2*x(5) + x(7)) - 2*l1*m1*m2*m4*x(2)^2*sin(x(1) - 2*x(5) + x(7)) - 4*l2*m1*m2*m3*x(4)^2*sin(x(3) - 2*x(5) + x(7)) - 4*l2*m1*m2*m4*x(4)^2*sin(x(3) - 2*x(5) + x(7)) - 2*l1*m1*m2*m3*x(2)^2*sin(x(1) + x(7)) - 4*l2*m1*m2*m3*x(4)^2*sin(x(3) + x(7)) - 8*l3*m1*m2*m3*x(6)^2*sin(x(5) + x(7)) + 2*l1*m1*m2*m3*x(2)^2*sin(x(1) - 2*x(3) + 2*x(5) - x(7)) + 2*l1*m1*m2*m4*x(2)^2*sin(x(1) - 2*x(3) + 2*x(5) - x(7)) - 2*l1*m1*m2*m3*x(2)^2*sin(2*x(3) - x(1) + x(7)) - 2*l1*m1*m2*m3*x(2)^2*sin(2*x(5) - x(1) + x(7)) - 4*l2*m1*m2*m3*x(4)^2*sin(2*x(5) - x(3) + x(7)))/(4*l4*m2^2*m3 + 2*l4*m2^2*m4 + 8*l4*m1*m2*m3 + 4*l4*m1*m2*m4 + 4*l4*m1*m3*m4 + 4*l4*m2*m3*m4 + l4*m2^2*m4*cos(2*x(1) - 2*x(3) - 2*x(5) + 2*x(7)) + l4*m2^2*m4*cos(2*x(1) - 2*x(3) + 2*x(5) - 2*x(7)) + 4*l4*m2^2*m3*cos(2*x(5)) - 2*l4*m2^2*m3*cos(2*x(3) - 2*x(1) + 2*x(5)) - 2*l4*m2^2*m3*cos(2*x(1) - 2*x(3) + 2*x(5)) - 4*l4*m2^2*m3*cos(2*x(1) - 2*x(3)) - 2*l4*m2^2*m4*cos(2*x(1) - 2*x(3)) - 2*l4*m2^2*m4*cos(2*x(5) - 2*x(7)) + 8*l4*m1*m2*m3*cos(2*x(5)) + 4*l4*m1*m3*m4*cos(2*x(7)) + 4*l4*m2*m3*m4*cos(2*x(7)) - 2*l4*m2*m3*m4*cos(2*x(3) - 2*x(1) + 2*x(7)) - 2*l4*m2*m3*m4*cos(2*x(1) - 2*x(3) + 2*x(7)) - 2*l4*m1*m3*m4*cos(2*x(5) - 2*x(3) + 2*x(7)) - 2*l4*m1*m3*m4*cos(2*x(3) - 2*x(5) + 2*x(7)) - 4*l4*m2*m3*m4*cos(2*x(1) - 2*x(3)) - 4*l4*m1*m3*m4*cos(2*x(3) - 2*x(5)) - 4*l4*m1*m2*m4*cos(2*x(5) - 2*x(7)));
        f = [f_1;f_2;f_3;f_4;f_5;f_6;f_7;f_8];
    end

X1 = l1*sin(x(:,1));
Y1 = l1*cos(x(:,1));
X2 = X1 + l2*sin(x(:,3));
Y2 = Y1 + l2*cos(x(:,3));
X3 = X2 + l3*sin(x(:,5));
Y3 = Y2 + l3*cos(x(:,5));
X4 = X3 + l4*sin(x(:,7));
Y4 = Y3 + l4*cos(x(:,7));
X0 = 0*X1;
Y0 = 0*Y1;
X = [X0 X1 X2 X3 X4];
Y = [Y0 Y1 Y2 Y3 Y4];

for n=1:10:50000

      P = plot(X(n,:),Y(n,:),'ro-');
      set(P,'LineWidth',1);
      set(P,'MarkerFaceColor','r');
      plot(X4(n),Y4(n),'b.');
      plot(X3(n),Y3(n),'g.');
      plot(X2(n),Y2(n),'c.');
      plot(X1(n),Y1(n),'m.');
      hold on
      pbaspect([1 1 1])
      xlim([-(l1+l2+l3+l4) (l1+l2+l3+l4)])
      ylim([-(l1+l2+l3+l4) (l1+l2+l3+l4)])
      pause(0.001)
      delete(P)

end






end