function LaufkatzeRK4 ()

% Parameter (global definiert, damit Wertaenderung nur an einer Stelle erforderlich):
global mK mL JL c lS Dt F0 a g
mK = 100 ; mL = 500 ; JL = 400 ; c = 200000 ; lS = 4 ; Dt = 1 ; F0 = 2000 ; a = 5 ; g = 9.81 ;

x0    = [0 ; 0 ; 0 ; 0] ;         % Anfangswerte
tspan = [0  10] ;                 % Zeitintervall

% Runge-Kutta 4. Ordnung: ---------------------------------
nrk4 = 5000 ;                   % Anzahl der Runge-Kutta-Schritte
[t xv] = rk4 (@RechteSeite , tspan , x0 , nrk4) ;

% Grafische Ausgabe von x und phi, zunaechst "Sortieren der Ergebnisse":
x     = xv(:,1) ;  % xv ist die Matrix der Ergebnisse,
phi   = xv(:,2) ;  % x, phi ... 
v     = xv(:,3) ;  % v und omega ... 
omega = xv(:,4) ;  % ... sind Vektoren 

Tges  = 0.5*mK*v.^2 + 0.5*mL*(v.^2+lS^2*omega.^2+2*lS*v.*omega.*cos(phi)) ...
       +0.5*JL*omega.^2-mL*g*lS*cos(phi) ;           % Kontrollfunktion "Gesamtenergie"   

subplot(3,1,1) ; plot (t , x)    , grid on , title  ('Bewegung der Laufkatze: x(t)') , axis([0 10 -25 10])
subplot(3,1,2) ; plot (t , phi)  , grid on , title  ('Pendeln der Last: \phi(t)')
subplot(3,1,3) ; plot (t , Tges) , grid on , title  ('Gesamtenergie: T_{ges}(t)')
disp(['Ergebnisse nach ' , num2str(length(t)-1) , ' Zeitschritten:'])  ; % "End"-Kontrolle
disp(['x    (t_end = ' , num2str(t(end)) , ') = ' , num2str(x  (end))])  ; 
disp(['phi  (t_end = ' , num2str(t(end)) , ') = ' , num2str(phi(end))])  ;
disp(['T_ges(t_end = ' , num2str(t(end)) , ') = ' , num2str(Tges(end))])  ;
disp(['(zum Vergleich: T_ges(t = ' , num2str(0) , ') = ' , num2str(Tges(1)) , ')'])  ;

% ============================================================
% Funktion, die die "rechte Seite" des Dgl.-Systems definiert: 
function r = RechteSeite (t , xv)

global mK mL JL c lS Dt F0 a g

x = xv(1) ; phi = xv(2) ; v = xv(3) ; omega = xv(4) ;

if (t < Dt) Ft = F0 ;  else   Ft = 0  ;  end
if (x < a ) ct = 0  ;  else   ct = c  ;  end

a11  = mK+mL ;
a12  = mL*lS*cos(phi) ;
a22  = mL*lS^2+JL ;
b1   = mL*lS*omega^2*sin(phi)-ct*(x-a)+Ft ;
b2   = -mL*g*lS*sin(phi) ;
detA = a11*a22-a12^2 ;

r = [v ; omega ; (b1*a22-b2*a12)/detA ; (b2*a11-b1*a12)/detA] ;