function [T,T1,fehler] = romberg(f,a,b,relerr)

% Funktion zur Romberg-Extrapolation mit vorgegebener
% Extrapolationsordnung
%
% Parameter:
% ----------
% a,b    : Integrationsintervallgrenzen
% f      : Integrand (als Zeiger in eigenen Funktionen, zB f1, f2, f3, f4)
% relerr : vorgegebene Fehlerordnung, es soll in einer Zeile zwei mal 
%          nacheinander abs(T_{ik}-T_{i-1,k})/T_{i,k}<=relerror gelten,
%          damit der Algorithmus abbricht und das Ergebnis T1 ausgibt. 
%
% Ausgaben: 
% ----------
% T1     : Wert des Integrals
% T      : gesamtes Tableau
% fehler : falls nach i=10 Zeilen immernoch nicht die gewuenschte
%          Genauigkeit erreicht wurde, ist es interessant zu wissen, welche
%          Genauigkeit erreicht wurde
%



% Startwert der Romberg Folge: 
% ----------------------------
h = (b-a);
hi(1) = h; 

% Indexshift T(0,0)->T(1,1)!! Trapezregel: 
% ----------------------------------------
T(1,1) = hi(1)*(f(a)+f(b))/2;

% Startwert fuer den Fehler:
% --------------------------
fehler = 1; 

% Index: 
% ------
i = 2; 
n = 1; 

while fehler>relerr
    
    % Romberg Folge: 
    % --------------
    hi(i) = hi(i-1)/2; % h wird in jedem Schritt halbiert      
    
    % Da h(i+1)=h(i)/2 gilt, folgt ohne weitere Funktionsauswertungen: 
    T(i,1) = T(i-1,1)/2;
    
    % Es muss in einem neuen Iterationsschritt nur jede zweite Auswertung
    % von f (und zwar die ungeraden) neu berechnet werden: 
    % (die Anzahl der Summanden verdoppelt sich in jedem Schritt)
    summe = 0;    
    for j=1:n
        summe = summe + f(a+(2*j-1)*hi(i));
    end
    T(i,1) = T(i,1) + hi(i)*summe;
    
    if i<=4
        j=i;
    else
        j=4; % hoechstens 4 Spalten berechnen... 
    end
    
    for k=2:j
        % rekursive Berechnung der T(i,k):
        T(i,k)=T(i,k-1)+(T(i,k-1)-T(i-1,k-1))/(2^(2*(k-1))-1);
        
        % Fehlerbestimmung:
        % (mind. 4 Zeilen und Spalten)
        if i>=4 & k>=4
            fehler1= abs(T(i,k)-T(i,k-1))/T(i,k);
            fehler2= abs(T(i,k-1)-T(i,k-2))/T(i,k-1);
            if fehler1<=relerr & fehler2<=relerr
                fehler = relerr;
            end
        end
    end
    
    % maximal 10 Spalten berechnen...: 
    % --------------------------------
    if i==10
        fehler=relerr;
    end
    
    % neue Indizes:
    i = i+1; 
    
    % da die Anzahl der Summanden sich in jedem Schritt verdoppelt:
    n=2*n;
end

T1 = T(i-1,k); % Wert des Integrals. Da i zuletzt einen erhoeht wurde, 
               % muss hier i-1 stehen.  