% Verschiedene Runge-Kutta-Verfahren fuer Stab-Pendel mit großen Ausschlaegen:

function ButchPendel (nsteps , verf)

if (nargin < 1) n = 100    ; Verfahren = 'RK4' ; 
else            n = nsteps ;
                if (nargin < 2) Verfahren = 'RK4' ;
                else            Verfahren = verf  ;
                end
end

global pl       g
       pl = 1 ; g = 9.81 ;                 % Pendellaenge, Erdbeschleunigung

t0   =  0 ;  phi0 = pi/2 ; omega0 = 0 ;    % Anfangswerte
tEnd = 10 ;                                % Ende des Integrationsintervalls

% Integrationsverfahren ---------------------------------------------------------- 
[t , xv] = butch (@BewDglPendel , [t0 ; tEnd] , [phi0 ; omega0] , n , Verfahren) ;
% Integrationsverfahren ----------------------------------------------------------

phi   = xv(:,1) ;                          % Sortieren der Ergebnisse:
omega = xv(:,2) ;                          % Spalte 1: phi, Spalte 2: omega

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

disp (['Integrationsverfahren:  ' , Verfahren]) ;
disp (['Anzahl der Zeitschritte n = ' , num2str(n)]) ;
disp (['phi_0     = ' , num2str(phi0)])   ; 
disp (['phi_end   = ' , num2str(phi(end),8)]) ;
disp (['omega_end = ' , num2str(omega(end),8)]) ; 
subplot (2,1,1) ; plot (t , phi)  , title ('\phi(t)') ;    
subplot (2,1,2) ; plot (t , Tges) , title ('T_{ges}(t)') ;

%%%%%%%%%% Function, die die Differenzialgleichungen definiert: %%%%%%%%%%%%%%%%%%%%%
function xvp = BewDglPendel (t , xv)
global pl g                    % ... sind in der Function DP45Pendel (oben) definiert.
phi    = xv(1) ;               % Input-Werte kommen in ...
omega  = xv(2) ;               % ... einem Vektor an.
phip   = omega ;               % Die beiden Differenzialgleichungen
omegap = -1.5*g/pl*sin(phi) ;  % werden ausgewertet, ...
xvp    = [phip ; omegap] ;     % ... die Ergebnisse werden im Vektor xvp abgeliefert.