function Autorad

% Parameter (global definiert, damit Wertaenderung nur an einer Stelle erforderlich):
global mK mR cK cR k a v0 b
mK = 1000 ; mR = 25 ; cK = 30000 ; cR = 90000 ; k = 750 ; a = 0.2 ; v0 = 50000/3600 ; b = 0.4 ;

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

% Integration des Anfangswertproblems:
maxstp  = 0.1
options = odeset('MaxStep' , maxstp) ; 
[t x]   = ode45 (@RechteSeite , tspan , x0 , options) ;

% Grafische Ausgabe, zunaechst "Sortieren der Ergebnisse":
xK  = [1 0 0 0] * x' ;  % x' ist die transponierte Matrix der Ergebnisse,
xR  = [0 1 0 0] * x' ;  
vK  = [0 0 1 0] * x' ;  
vR  = [0 0 0 1] * x' ;  

tlen = length(t) ;
xKEnd = xK(tlen)
xREnd = xR(tlen)
for ii = 1:tlen
    if (t(ii) < b/v0)   u(ii) = a/b*v0*t(ii) ;  else  u(ii) = a  ;  end
end

subplot(3,1,1) ; plot (t , xK ,'-r' , t , xR ,'-b' , t , u ,'-g') , 
                 grid on , title  ('x-t-Diagramm:') , legend ('xK','xR','u')
subplot(3,1,2) ; plot (t , vK ,'-r' , t , vR, '-b') , 
                 grid on , title  ('v-t-Diagramm:') , legend ('vK','vR')
subplot(3,1,3) ; plot (xK , vK) , grid on , title  ('vK-xK-Diagramm:')

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

global mK mR cK cR k a v0 b

xK  = x(1) ;  
xR  = x(2) ;  
vK  = x(3) ;  
vR  = x(4) ;  

if (t < b/v0)   u = a/b*v0*t ;  else  u = a  ;  end

xKp  = vK ;  
xRp  = vR ;  
vKp  = (-cK*(xK-xR)-k*(vK-vR))/(mK/4) ;  
vRp  = (cK*(xK-xR)-cR*(xR-u)+k*(vK-vR))/mR ;  

f = [xKp ; xRp ; vKp ; vRp] ;
