% Runge-Kutta-Verfahren 4. Ordnung fuer Stab-Pendel mit großen Ausschlaegen:

phip   = inline ('omega','omega') ;
omegap = inline ('-faktor*sin(phi)','t','phi','omega','faktor') ;

pl     = 1 ; g = 9.81 ;                % Pendellaenge, Erdbeschleunigung
faktor = 1.5*g/pl ;

t0 = 0 ; phi0 = pi/2 ; omega0 = 0 ;    % Anfangswerte
tEnd = 10   ;                          % Ende des Integrationsintervalls
n    = 1000 ;                          % Anzahl der Integrationsschritte
dt   = (tEnd - t0) / n ;

t = t0 : dt : tEnd    ;
phi   = zeros (n+1,1) ;
omega = zeros (n+1,1) ;

% Runge-Kutta-Algorithmus -------------------------------------------------- 
phi(1) = phi0 ; omega(1) = omega0 ;

for i=1:n
    k1phi = phip   (omega(i)) ;
    k1om  = omegap (t(i),phi(i),omega(i),faktor) ;
    k2phi = phip   (omega(i)+k1om*dt/2) ;
    k2om  = omegap (t(i)+dt/2,phi(i)+k1phi*dt/2,omega(i)+k1om*dt/2,faktor) ;
    k3phi = phip   (omega(i)+k2om*dt/2) ;
    k3om  = omegap (t(i)+dt/2,phi(i)+k2phi*dt/2,omega(i)+k2om*dt/2,faktor) ;
    k4phi = phip   (omega(i)+k3om*dt) ;
    k4om  = omegap (t(i+1),phi(i)+k3phi*dt,omega(i)+k3om*dt,faktor) ;
    phi  (i+1) = phi  (i) + (k1phi + 2*k2phi + 2*k3phi + k4phi) * dt/6 ;
    omega(i+1) = omega(i) + (k1om  + 2*k2om  + 2*k3om  + k4om)  * dt/6  ;
end
% Runge-Kutta-Algorithmus --------------------------------------------------

% Erweiterte Auswertung, Verifizieren der Ergebnisse:
Tges = - g/2*cos(phi) + pl*omega.^2/6 ;   % Kontrollfunktion "Gesamtenergie"

[phimin , imin] = min(phi) ; 
disp (['Anzahl der Zeitschritte n = ' , num2str(n)]) ;
disp (['phi_0     =  ' , num2str(phi0)])   ; 
disp (['phi_min   = ' , num2str(phimin) , ' bei t = ' , num2str(t(imin))]) ;
disp (['omega_max =  ' , num2str(max(abs(omega)))])   ; 
subplot (2,1,1) ; plot (t , phi)  , title ('\phi(t)') ;    
subplot (2,1,2) ; plot (t , Tges) , title ('T_{ges}(t)') ;      