% Stab-Pendel mit großen Ausschlaegen, Loesung mit Matlab-Standard-Solver ode45:

function DP45oPendel (epstol , mstep) 

if (nargin < 1) reltol = 1e-3 ; maxstep = 1 ; 
else            reltol = epstol ;
                if (nargin < 2) maxstep = 1 ; 
                else            maxstep = mstep ;
                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

% Dormand-Prince-Algorithmus ----------------------------------------------- 
options  = odeset ('RelTol' , reltol , 'MaxStep' , maxstep) ;
[t , xv] = ode45 (@BewDglPendel , [t0 ; tEnd] , [phi0 ; omega0] , options) ;
% Dormand-Prince-Algorithmus ----------------------------------------------- 

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"

[phimin , imin] = min(phi) ; 
disp (['Anzahl der Zeitschritte n = ' , num2str(length(t))]) ;
disp (['eps_rel   = ' , num2str(reltol)])   ; 
disp (['MaxStep   = ' , num2str(maxstep)])   ; 
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 RKPendel4 (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.