% Die Funktion rk4 loest ein Anfangswertproblem mit den Runge-Kutta-Formeln
% 4. Grades. Es wird mit der konstanten Schrittweite (tspan(2)-tspan(1))/n
% gerechnet, so dass auf t ein aequidistantes Feld mit n+1 Werten abgeliefert wird.
% Die Anzahl der zu loesenden Differenzialgleichungen 1. Ordnung wird durch die 
% Laenge des Vektors der Anfangsbedingungen x0 festgelegt, die Differenzialgleichungen
% selbst muessen in einer Funktion definiert sein, deren Handle funh an rk4
% uebergeben werden muss.
%
% Auf der Matrix x (n+1 Zeilen) werden die Ergebnisse der numerischen
% Integration fuer alle in t verzeichneten Zeitpunkte abgeliefert.
% Die output-Werte t und x und die ersten 3 input-Parameter funh, tspan, x0
% entsprechen genau denen der MATLAB-ode-Funktionen, n ist die feste Anzahl
% der auszufuehrenden Integrationsschritte.

function [t , x] = rk4 (funh , tspan , x0 , n)

tanf = tspan(1) ;
tend = tspan(2) ;
dt   = (tend - tanf) / n ;
ndeq = length (x0) ;

t = tanf : dt : tend ;
x = zeros (n+1,ndeq) ;

% Runge-Kutta-Algorithmus -------------------------------------------------- 
for i=1:ndeq 
    x(1,i) = x0(i) ;
end

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