L1 =1;L2=1;m1=1;m2=2;g=9.8;
tspan = [0:1/30:20];
dpend = @(t,x) [x(2); (-m2*L1*(x(2))^2*sin(x(1)-x(3))*cos(x(1)-x(3)) + ...
    g*m2*sin(x(3))*cos(x(1)-x(3)) - m2*L2*(x(4))^2*sin(x(1)-x(3))-(m1+m2)...
    *g*sin(x(1)))/(L1*(m1+m2) - m2*L1*(cos(x(1)-x(3)))^2); x(4); (m2*L2*...
    (x(4))^2*sin(x(1)-x(3))*cos(x(1)-x(3))+g*sin(x(1))*cos(x(1)-x(3))*...
    (m1+m2)+L1*(x(2))^2*sin(x(1)-x(3))*(m1+m2)-g*sin(x(3))*(m1+m2))/(L2*...
    (m1+m2)-m2*L2*(cos(x(1)-x(3)))^2)];


initial2 = [4.43-0.001;0;2.75+0.001;0];
[t,x] = ode45(dpend,tspan,initial2);

figure
% Position of mass 1 with respect to the origin
x1 = L1*sin(x(:,1)); y1 = -L1*cos(x(:,1));
% Position of mass 2 with respect to the origin
x2 = x1+L2*sin(x(:,3)); y2 = y1-L2*cos(x(:,3));

itmax= 280;
f = 0;
for iter = 1:itmax
    clf;
    subplot(1,2,1)
    line([0 x1(iter)],[0 y1(iter)]); hold on;
    subplot(1,2,1)
    plot(x1(iter),y1(iter),'ro');
    subplot(1,2,1)
    line([x1(iter) x2(iter)],[y1(iter) y2(iter)]);
    subplot(1,2,1)
    plot(x2(iter),y2(iter),'ro');
    title(sprintf('Motion at time: %.1f seconds',t(iter))); xlabel('x'); ylabel('y');
    axis equal; axis([-2,2,-2,2]); hold off; pause(1/30); drawnow;
    
    subplot(1,2,2)
    plot(x2(1:iter),y2(1:iter),'k--');
    title(sprintf('Trace of double pendulum @ time: %.1f seconds',t(iter))); xlabel('x'); ylabel('y');
    axis equal; axis([-2,2,-2,2]); hold off; pause(1/30); drawnow;
    
    f=f+1;
    saveas(gcf, ['movie2',sprintf('%03d',f),'.png'])
    pause(0.01)
    
end