% Differenzialgleichungen:

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

pl     = 1 ;                               % Pendellaenge in m
g      = 9.81 ;                            % Erdbeschleunigung in m/s^2
faktor = 1.5*g/pl ;

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

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

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

for i=1:n
    phipi = phip   (omega(i)) ;
    ompi  = omegap (t(i),phi(i),omega(i),faktor) ;
    % * Prädiktorschritt nach Euler-Cauchy: ********************************
    prphi = phi  (i) + dt * phipi ;
    prom  = omega(i) + dt * ompi  ;
    for j=1:3   % * Drei Korrektorschritte: ********************************
        prphi = phi  (i) + (phipi + phip(prom))*dt/2 ;
        prom  = omega(i) + (ompi  + omegap(t(i+1),prphi,prom,faktor))*dt/2 ;
    end
    phi  (i+1) = prphi ;           % ... die endgueltigen Werte des Schritts
    omega(i+1) = prom  ;
end
% Heun-Algorithmus ---------------------------------------------------------

phiMin=min(phi) , omegaMax=max(abs(omega)) % Extremwerte ---> Command Window
plot (t,phi)                               % Graphik-Ausgabe phi(t)
