% Die Funktion butch loest ein Anfangswertproblem mit einem Runge-Kutta-Verfahren
% 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 xv0 festgelegt, die Differenzialgleichungen
% selbst muessen in einer Funktion definiert sein, deren Handle funh an butch
% uebergeben werden muss.
%
% Auf der Matrix xv (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, xv0
% entsprechen genau denen der MATLAB-ode-Funktionen, n ist die feste Anzahl
% der auszufuehrenden Integrationsschritten but ist das Runge-Kutta-Verfahren,
% das fuer die Integration verwendet werden soll. Folgende Verfahren sind
% verfuegbar:
%
% but = 'EC'        --->  Verfahren von Euler-Cauchy
% but = 'Heun'      --->  Verfahren von Heun mit Euler-Cauchy als Praediktor
%                         (keine Iteration)
% but = 'RK4'       --->  Klassisches Runge-Kutta-Verfahren 4. Ordnung
% but = 'Fehlberg4' --->  Verfahren von Runge-Kutta-Fehlberg 4. Ordnung
% but = 'Fehlberg5' --->  Verfahren von Runge-Kutta-Fehlberg 5. Ordnung
% but = 'DP4'       --->  Verfahren von Dormand-Prince 4. Ordnung
% but = 'DP5'       --->  Verfahren von Dormand-Prince 5. Ordnung

function [t , xv] = butch (funh , tspan , xv0 , n , but)

ndeq = length (xv0) ;

if (strcmp(but, 'Heun'))          [s , c , b , A] = Heun ;
elseif (strcmp(but , 'EC'))       [s , c , b , A] = EC   ;
elseif (strcmp(but , 'RK4'))      [s , c , b , A] = RK4  ;
elseif (strcmp(but, 'Fehlberg4')) [s , c , b , A] = Fehlberg (4) ;
elseif (strcmp(but, 'Fehlberg5')) [s , c , b , A] = Fehlberg (5) ;
elseif (strcmp(but, 'DP4'))       [s , c , b , A] = DP (4) ;
elseif (strcmp(but, 'DP5'))       [s , c , b , A] = DP (5) ;
else   error('Fehler: Unbekanntes Integrationsverfahren')
end

k = zeros (s , ndeq) ;
tanf = tspan(1) ;
tend = tspan(2) ;
dt   = (tend - tanf) / n ;

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

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

for i=1:n
    ti = t(i) ;
    for j=1:s
        tij  = ti + dt*c(j) ;
        xvij = xv(i,:) ;
        for l=1:j-1
            xvij = xvij + dt*A(j,l)*k(l,:) ;
        end
        k(j,:) = feval (funh , tij , xvij) ;
    end
    xv(i+1,:) = xv(i,:) ;
    for j=1:s
        xv(i+1,:) = xv(i+1,:) + dt*b(j)*k(j,:) ;
    end
end
% Runge-Kutta-Algorithmus
% --------------------------------------------------

function [s , c , b , A] = EC     
    s = 1 ;
    c = [ 0 ] ;
    b = [ 1 ] ;
    A = [ 0 ] ;

function [s , c , b , A] = Heun     
    s = 2 ;
    c = [ 0   1 ]   ;
    b = [ 1/2 1/2 ] ;
    A = [ 0   0   ; 
          1   0   ] ;

function [s , c , b , A] = RK4     
    s = 4 ;
    c = [ 0   1/2 1/2 1 ]   ;
    b = [ 1/6 1/3 1/3 1/6 ] ;
    A = [ 0   0   0   0     ; 
          1/2 0   0   0     ;
          0   1/2 0   0     ;
          0   0   1   0   ] ;
      
function [s , c , b , A] = Fehlberg (Ordnung)     
    s = 6 ;
    c = [     0         1/4        3/8      12/13      1  1/2] ;
    if (Ordnung == 4)
        b = [ 25/216  0  1408/2565   2197/4104  -1/5   0   ] ;
    else
        b = [ 16/135  0  6656/12825 28561/56430 -9/50 2/55 ] ;
    end
    A = [     0          0          0         0        0   0   ;
             1/4         0          0         0        0   0   ;
             3/32       9/32        0         0        0   0   ;
          1932/2197 -7200/2197  7296/2197     0        0   0   ;
           439/216      -8      3680/513  -845/4104    0   0   ;
            -8/27        2     -3544/2565 1859/4104 -11/40 0 ] ;

function [s , c , b , A] = DP (Ordnung)     
    s = 7 ;
    c = [     0         1/5        3/10      4/5  8/9  1  1] ;
    if (Ordnung == 4)
        b = [ 5179/57600 0 7571/16695 393/640 -92097/339200 187/2100 1/40 ] ;
    else
        b = [   35/384   0  500/1113  125/192  -2187/6784    11/84    0   ] ;
    end
    A = [     0           0          0         0         0        0   0   ;
             1/5          0          0         0         0        0   0   ;
             3/40        9/40        0         0         0        0   0   ;
            44/45      -56/15      32/9        0         0        0   0   ;
         19372/6561 -25360/2187 64448/6561 -212/729      0        0   0   ;
          9017/3168   -355/33   46732/5247   49/176 -5103/18656   0   0   ;
            35/384        0       500/1113  125/192 -2187/6784  11/84 0 ] ;
