View difference between Paste ID: adDfzz96 and 87cuwG6S

 function quadruplependulumsim
-close all
+    t = 0:0.001:60;
-t=0:0.001:60;
+    g = 9.81;
-g=9.81;
+    m1 = 1;
-m1=1;
+    m2 = 1;
-m2=1;
+    m3 = 1;
-m3=1;
+    m4 = 1;
-m4=1;
+    l1 = 1.6;
-l1=1.6;
+    l2 = 1.4;
-l2=1.4;
+    l3 = 1.2;
-l3=1.2;
+    l4 = 1;
-l4=1;
+
-theta1_initial = 0.2;
+    % initial solution: [theta1, theta1dot, theta2, theta2dot, theta3, theta3dot, theta4, theta4dot]
-theta1dot_initial = 0;
+    x0 = [0.2, 0, 0.4, 0, 0.6, 0, 0.8, 0];
-theta2_initial = 0.4;
+    [t, x] = ode45(@rhs, t, x0);
-theta2dot_initial = 0;
+
-theta3_initial = 0.6;
+    figure(1)
-theta3dot_initial = 0;
+    plot(t, x(:,[1 3 5 7]))
-theta4_initial = 0.8;
+    xlabel('t (sec)')
-theta4dot_initial = 0;
+    ylabel('Angle (rad)')
-[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));
+    l = [l1, l2, l3, l4];
-xlabel('t (s)'),ylabel('Angle (rad)');
+    X = cumsum(bsxfun(@times, sin(x(:,[1 3 5 7])), l), 2);
-%pause
+    Y = cumsum(bsxfun(@times, cos(x(:,[1 3 5 7])), l), 2);
-    function f = rhs(t,x)
+
     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{:});
-        f = [f_1;f_2;f_3;f_4;f_5;f_6;f_7;f_8];
+        else
             hLines(i) = line(NaN, NaN, 'Color',clr(i), opts{:});
         end
-X1 = l1*sin(x(:,1));
+
-Y1 = l1*cos(x(:,1));
+    hLine = line(NaN, NaN, ...
-X2 = X1 + l2*sin(x(:,3));
+        'LineWidth',1, 'Marker','o', 'Color','r', 'MarkerFaceColor','r');
-Y2 = Y1 + l2*cos(x(:,3));
+    axis equal
-X3 = X2 + l3*sin(x(:,5));
+    axis([-1 1 -1 1] * sum(l))
-Y3 = Y2 + l3*cos(x(:,5));
+
-X4 = X3 + l4*sin(x(:,7));
+    for n=1:10:size(X,1)
-Y4 = Y3 + l4*cos(x(:,7));
+        for i=1:4
-X0 = 0*X1;
+            if isHG2
-Y0 = 0*Y1;
+                addpoints(hLines(i), X(n,i), Y(n,i));
-X = [X0 X1 X2 X3 X4];
+            else
-Y = [Y0 Y1 Y2 Y3 Y4];
+                set(hLines(i), 'XData',X(1:10:n,i), 'YData',Y(1:10:n,i));
             end
-for n=1:10:50000
+        end
         set(hLine, 'XData',[0 X(n,:)], 'YData',[0 Y(n,:)]);
-      P = plot(X(n,:),Y(n,:),'ro-');
+        drawnow
-      set(P,'LineWidth',1);
+
-      set(P,'MarkerFaceColor','r');
+
-      plot(X4(n),Y4(n),'b.');
+    function f = rhs(t, x)
-      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)
+        f = [f_1; f_2; f_3; f_4; f_5; f_6; f_7; f_8];
     end
-end
+
 end