quadruple pendulum motion

function quadruplependulumsim
    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;

    % initial solution: [theta1, theta1dot, theta2, theta2dot, theta3, theta3dot, theta4, theta4dot]
    x0 = [0.2, 0, 0.4, 0, 0.6, 0, 0.8, 0];
    [t, x] = ode45(@rhs, t, x0);

    figure(1)
    plot(t, x(:,[1 3 5 7]))
    xlabel('t (sec)')
    ylabel('Angle (rad)')

    l = [l1, l2, l3, l4];
    X = cumsum(bsxfun(@times, sin(x(:,[1 3 5 7])), l), 2);
    Y = cumsum(bsxfun(@times, cos(x(:,[1 3 5 7])), l), 2);

    isHG2 = ~verLessThan('matlab', '8.4');

    figure(2)
    clr = 'mcgb';
    opts = {'LineStyle','none', 'Marker','.', 'MarkerSize',1}; % 6
    for i=4:-1:1
        if isHG2
            hLines(i) = animatedline('Color',clr(i), opts{:});
        else
            hLines(i) = line(NaN, NaN, 'Color',clr(i), opts{:});
        end
    end
    hLine = line(NaN, NaN, ...
        'LineWidth',1, 'Marker','o', 'Color','r', 'MarkerFaceColor','r');
    axis equal
    axis([-1 1 -1 1] * sum(l))

    for n=1:10:size(X,1)
        for i=1:4
            if isHG2
                addpoints(hLines(i), X(n,i), Y(n,i));
            else
                set(hLines(i), 'XData',X(1:10:n,i), 'YData',Y(1:10:n,i));
            end
        end
        set(hLine, 'XData',[0 X(n,:)], 'YData',[0 Y(n,:)]);
        drawnow
    end

    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

end