function DoppelPendelAnimation (phi1Anf , phi2Anf)

if nargin == 2                          % Start mit vorgegebenen Anfangswerten?
   x0 = [phi1Anf ; phi2Anf ; 0 ; 0] ;   % Anfangswerte  [phi1 ; phi2 ; omega1 ; omega2]
else
   x0 = [pi/2 ; 0 ; 0 ; 0] ;            % Default-Anfangswerte
end

% Parameter (global definiert, damit Wertaenderung nur an einer Stelle erforderlich):
global mdm  J1  J2  s1  s2  gdl 
mdm = 1 ; J1 = 1./12 ; J2 = 1./12 ; s1 = .5 ; s2  = .5 ; l1 = 1 ; l2 = 1 ; gdl = 9.81/l1 ; 

tspan = [0  10] ;                       % Zeitintervall

% Integration des Anfangswertproblems:
clf
options = odeset('Mass' , @Massenmatrix , 'MaxStep' , 0.001 , ...
                 'OutputFcn' , @odeplot , 'OutputSel' , [2]) ;
% ... zeichnet die Funktion phi2(t) waehrend der ode45-Berechnung in das Graphik-Fenster             
[t x]   = ode45 (@RechteSeite , tspan , x0 , options) ;

% Animation der Bewegung, zunaechst "Sortieren der Ergebnisse":
phi1i = x(:,1) ;    % x ist die Matrix der Ergebnisse,
phi2i = x(:,2) ; ;  % phi1i und phi2i sind Vektoren 

clf
tt = tspan(1) ;
xx = [0  l1*sin(phi1i(1))  l1*sin(phi1i(1))+l2*sin(phi2i(1))] ;
yy = [0 -l1*cos(phi1i(1)) -l1*cos(phi1i(1))-l2*cos(phi2i(1))] ;
plot ([0 l1/8 -l1/8 0],[0 l1/6 l1/6 0] , 'r-') ;
hold on
p=plot (xx,yy,'o-','EraseMode','xor') ;
axis (1.1*[-(l1+l2) (l1+l2) -(l1+l2) (l1+l2)]) ;
zeit = title (sprintf('t = %8.2f s' , tt)) ;
set (gca , 'userdata' , zeit)

deltatt  = 1/25     ;    % 25 Bilder pro Sekunde
zeitlupe = 1        ;    % Verzögerungsfaktor
cstart   = cputime  ;
Zeitschritte = length (t)

for ii=1:Zeitschritte    % Schleife über alle berechneten Zeitschritte
   while (t(ii) >= tt)   % Ausgabe jeweils nur nach deltatt 
       xx = [0  l1*sin(phi1i(ii))  l1*sin(phi1i(ii))+l2*sin(phi2i(ii))] ;
       yy = [0 -l1*cos(phi1i(ii)) -l1*cos(phi1i(ii))-l2*cos(phi2i(ii))] ;
       set(p , 'XData' , xx , 'YData' , yy)
       while (cputime-cstart < tt*zeitlupe)  % Warteschleife
       end 
       zeit = get (gca , 'userdata');
       set(zeit , 'string' , sprintf('t = %8.2f s' , tt))
       drawnow
       tt = tt + deltatt ;
   end
end

% ============================================================
% Funktion, die die "Massenmatrix" des Dgl.-Systems definiert: 
function M = Massenmatrix (t , xv)

global mdm  J1  J2  s1  s2  gdl 

phi1   = xv(1) ;
phi2   = xv(2) ;
omega1 = xv(3) ;
omega2 = xv(4) ;

c = cos(phi1-phi2);
M = [ 1  0      0             0      ; 
      0  1      0             0      ; 
      0  0  s1^2+J1+mdm    mdm*s2*c  ;
      0  0    mdm*s2*c   mdm*s2^2+J2 ] ;

% ============================================================
% Funktion, die die "rechte Seite" des Dgl.-Systems definiert: 
function xvp = RechteSeite (t , xv)

global mdm  J1  J2  s1  s2  gdl 

phi1   = xv(1) ;
phi2   = xv(2) ;
omega1 = xv(3) ;
omega2 = xv(4) ;

s = sin(phi1-phi2);

phi1p   = omega1 ;
phi2p   = omega2 ;
omega1p = -mdm*s2*omega2^2*s - (s1+mdm)*gdl*sin(phi1) ;
omega2p =  mdm*s2*omega1^2*s -   mdm*s2*gdl*sin(phi2) ;

xvp = [phi1p ; phi2p ; omega1p ; omega2p] ; 