
function [R,error] = romberg(F,a,b,n,tol)

% function [R,error] = romberg(F,a,b,n,tol)
% Berechnung des Integrals R = int_a^b F
% einer vektorwertigen Funktion [f_1; f_1; ...] 
% absolute Toleranz tol, mindestens n Stuetzstellen
% error: geschaetzter Fehler 
% MAX_STEPS = 8: maximale Anzahl der Extrapolationen

MAX_STEPS = 8;
h = (b-a)/(n-1); 

% erste Trapezsumme, Spaltenvektor der Dimension von F
x = [a:h:b+h/2];
Fx = feval(F,x); 
Fx(:,1) = Fx(:,1)/2; 
Fx(:,n) = Fx(:,n)/2;
S = sum(Fx')'*h; R = S;

% Extrapolation
for j=1:MAX_STEPS;
   % j-te Trapezsumme
   x = [a+h/2:h:b];
   Fx = feval(F,x);
   % letze Tableau-Zeile, ergaenzt um Trapezsumme
   S = [(S(:,1)+h*sum(Fx')')/2 S];
   % Ueberscheiben mit neuer Tableau-Zeile
   for k=1:j;
      S(:,k+1) = (4^k*S(:,k) - S(:,k+1))/(4^k-1);
   end;
   h = h/2; R = S(:,j+1);
   % Fehlerschaetzung 
   error = norm(abs(S(:,j)-S(:,j+1)),inf),
   if error < tol;
      return;
   end;
end;

display('MAX_STEPS extrapolations exceeded');


