function DoublePendulum()
t  = 0:0.01:10;
Y0 = [pi/2;pi/2;-2;0];
%Y0 = [2*pi/2,2*pi/2,-1,0];

function dy = DoubPend(y)
  L = 1;  g = 9.81;
  ddy = [2, cos(y(1)-y(2));
         cos(y(1)-y(2)),1]\(sin(y(1)-y(2))*[-y(4)^2; y(3)^2]...
			   -g/L*[2*sin(y(1));sin(y(2))]);
  dy = [y(3);y(4);ddy];
end%function

vopt = odeset('RelTol', 1e-8, 'AbsTol', 1e-8,'NormControl', 'on');
[t,Y] = ode45(@(t,y)DoubPend(y),t,Y0,vopt);

figure(1);clf
plot(t,Y(:,1)*180/pi,t,Y(:,2)*180/pi);
grid on; xlabel('time t'); ylabel('angles [deg]');
legend('first angle \phi','second angle \theta');
drawnow()

figure(2)
L = 1;
for ii = 1:length(t)
  data = Y(ii,1:2);
  plot([0 L*sin(data(1)) L*sin(data(1))+L*sin(data(2))],...
       [0 -L*cos(data(1)) -L*cos(data(1))-L*cos(data(2))],'o-')
  axis(2.1*[-L,L,-L,L])
  axis 'square'
  pause(0.01)
end%for
end
