function pendel ()

global pl g
          g = 9.81 ;

plbereich = [0.01 1] ;
nsteps    = 30 ;
plstep    = (plbereich(2) - plbereich(1))/nsteps ;
plen = zeros (nsteps+1 , 1) ;
tend = zeros (nsteps+1 , 1) ;


phi0  = [3*pi/4 ; 0] ;      % Anfangswerte
tspan = [0 1000]     ;      % Zeitintervall (groß genug)

% Optionen für Integration des Anfangswertproblems:
options = odeset ('MaxStep' , 0.01 , 'Events' , @Ereignis) ;

pl = plbereich(1) ;
for i=1:nsteps+1
   [t xv te xve ereig]  = ode45 (@BewDglPendel , tspan , phi0 , options) ;
   plen(i) = pl          ;
   tend(i) = t(end) * 2   ;     % Zeit für eine volle Schwingung
   pl      = pl + plstep ;
end

plot (plen , tend) , grid on , title ('pl-T-Diagramm:')

[plen tend]        %  ... schreibt Wertetabelle in das Command Window

% ===================================================================
function [value,isterminal,direction] = Ereignis (t,xvt)
global pl g

omega = xvt(2) ;
value(1)      = omega ; % Ereignis 1: "Winkelgeschwindigkeit Null"
isterminal(1) = 1 ;     % Abbruch, ...
direction(1)  = 0 ;     % ... in jedem Fall

% ===================================================================
% Funktion, die die "rechte Seite" des Dgl.-Systems definiert:

function xvp = BewDglPendel (t , xv)
global pl g                    % ... sind in 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.