% Stabpendel mit Runge-Kutta, Differenzialgleichungen in Function File
% BewDglPendel, angepasst an die Form, die die MATLAB-ode-Functions erwarten

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

t0 = 0 ;  phi0 = pi/2 ; omega0 = 0 ;       % Anfangswerte

tEnd = 10 ;                                % Ende des Integrationsintervalls
     
n    = 200 ;                               % Anzahl der Integrationsschritte
dt   = (tEnd - t0) / n ;

t = t0 : dt : tEnd    ;
x = zeros (n+1,2) ;

% Runge-Kutta-Algorithmus -------------------------------------------------- 
x(1,1) = phi0   ;
x(1,2) = omega0 ;

for i=1:n
    k1 = BewDglPendel (t(i),x(i,:)) ;
    k2 = BewDglPendel (t(i)+dt/2,x(i,:)+k1'*dt/2) ;
    k3 = BewDglPendel (t(i)+dt/2,x(i,:)+k2'*dt/2) ;
    k4 = BewDglPendel (t(i+1),x(i,:)+k3'*dt) ;
    x(i+1,:) = x (i,:) + (k1' + 2*k2' + 2*k3' + k4') * dt/6 ;
end
% Runge-Kutta-Algorithmus --------------------------------------------------

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