function DoppelPendelKontr (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 ; gdl = 9.81 ; 

tspan = [0 ; 10] ;                       % Zeitintervall

% Integration des Anfangswertproblems:
options = odeset('Mass' , @Massenmatrix , 'MaxStep' , 0.1) ; 
[t xv]  = ode45 (@RechteSeite , tspan , x0 , options) ;
[n nn] = size(t) ;                       % Anzahl der Zeitschritte

% Runge-Kutta-Algorithmus: -------------------------------------------------- 
[trk xvrk] = rk4 (@DopendelDgln , tspan , x0 , n) ;

% Grafische Ausgabe von phi1 und phi2, zunaechst "Sortieren der Ergebnisse":
phi1 = xv(:,1) ; 
phi2 = xv(:,2) ;                         % phi1 und phi2 sind Vektoren 
phi1rk = xvrk(:,1) ; 
phi2rk = xvrk(:,2) ;                     % phi1i und phi2i sind Vektoren 

subplot(3,1,1) ; plot (t , phi1 , 'r' , trk , phi1rk , 'k') , grid on , title  ('\phi_1(t)')
subplot(3,1,2) ; plot (t , phi2 , 'r' , trk , phi2rk , 'k') , grid on , title  ('\phi_2(t)')

Tq = zeros (n,1) ;
Tqrk = zeros (n,1) ;
om1   = xv(:,3) ; 
om2   = xv(:,4) ;                   
om1rk = xvrk(:,3) ; 
om2rk = xvrk(:,4) ;                   
for i=1:n
    Tq(i)   = 0.5*(s1^2+J1+mdm)/gdl*om1(i)^2 + 0.5*(mdm*s2^2 + J2)/gdl*om2(i)^2 + (mdm*s2*cos(phi1(i)-phi2(i)))/gdl*om1(i)*om2(i) - (s2+mdm)*cos(phi1(i)) - mdm*s2*cos(phi2(i)) ; 
    Tqrk(i) = 0.5*(s1^2+J1+mdm)/gdl*om1rk(i)^2 + 0.5*(mdm*s2^2 + J2)/gdl*om2rk(i)^2 + (mdm*s2*cos(phi1rk(i)-phi2rk(i)))/gdl*om1rk(i)*om2rk(i) - (s2+mdm)*cos(phi1rk(i)) - mdm*s2*cos(phi2rk(i)) ; 
end
subplot(3,1,3) ; plot (t , Tq , 'r' , t , Tqrk , 'k') , grid on , title  ('T_{ges}(t)')

disp(['Ergebnisse nach ' , num2str(length(t)-1) , ' Zeitschritten:'])  ; % "End"-Kontrolle
disp(['ode45: phi1 (t_end = ' , num2str(t(end)) , ') = ' , num2str(phi1(end))])  ; 
disp(['ode45: phi2 (t_end = ' , num2str(t(end)) , ') = ' , num2str(phi2(end))])  ;
disp(['ode45: |T|_max = ' , num2str(max(abs(Tq))) , '  |T|_min = ' , num2str(min(abs(Tq)))])  ;
disp(['rk4:   phi1 (t_end = ' , num2str(t(end)) , ') = ' , num2str(phi1rk(end))])  ; 
disp(['rk4:   phi2 (t_end = ' , num2str(t(end)) , ') = ' , num2str(phi2rk(end))])  ;
disp(['rk4:   |T|_max = ' , num2str(max(abs(Tqrk))) , '  |T|_min = ' , num2str(min(abs(Tqrk)))])  ;

% ============================================================
% 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)   % ... fuer ode45

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] ;

% ============================================================
% Funktion, die die Differenzialgleichungen definiert: 
function xvp = DopendelDgln (t , xv)   % ... fuer rk4

global mdm  J1  J2  s1  s2  gdl 

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

s = sin(phi1-phi2) ;

a11 = s1^2 + J1 + mdm ;
a12 = mdm * s2 * cos(phi1-phi2) ;
a22 = mdm * s2^2 + J2 ;
b1  = -mdm*s2*omega2^2*s - (s1+mdm)*gdl*sin(phi1) ;
b2  =  mdm*s2*omega1^2*s -   mdm*s2*gdl*sin(phi2) ;

detA = a11 * a22 - a12^2 ;

phi1p   = omega1 ;
phi2p   = omega2 ;
omega1p = (b1 * a22 - b2 * a12)/detA ; 
omega2p = (b2 * a11 - b1 * a12)/detA ; 

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