function [result, T] = romberg(intfunc, quadratur, a, b, N, m)
%
%  Berechnung des Rombergschemas fuer den Integranden integrand
%  mit der Quadraturformel quadratur.
%  (quadratur koennte ein allgemeiner Operator sein)
%
%  intfunc    Funktion der Form "function y=integrand(x)",
%              zB. 'func_a'
%  quadratur  Quadratur der Form
%             "function I=quadratur(integrand, a, b, N)",
%             zB. 'trapez'
%  a          untere Intervallgrenze
%  b          obere  Intervallgrenze
%  N+1        Anzahl Stuetzstellen (aequdistent) auf der groebsten Stufe
%  m          Anzahl Stufen
%  T          Romberg-Schema
%  result     der 'beste' Werte fuer das Integral
%
%  Aufruf der Funktion z.B.:
%    [result, T] = romberg('func_a', 'trapez', 0, 3, N, m);
%    (mit sinnvollen Werten fuer N und m)
%
%  Hinweis:
%  Das Programm definiert eine leere Matrix T fuer das Romberg-Schema.
%  Zunaechst wird die erste Spalte gefuellt:
%  Das Element T(1,1) ist das Resultat der Quadratur mit  N Teilintervallen,
%  das Element T(2,1) ist das Resultat der Quadratur mit 2N Teilintervallen,
%  das Element T(3,1) ist das Resultat der Quadratur mit 4N Teilintervallen,
%  usw. bis zum Element T(m+1,1)
%  Die Aufgabe besteht darin, an der bezeichneten Stelle den Code fuer
%  die Romberg-Extrapolation einzufuegen.
%  Dabei wird man Spaltenweise vorgehen (eine for-Schleife fuer die Spalten
%  2 bis m+1) und eine innere for-Schleife fuer die Zeilen von T.

% Romberg-Schema initialisieren
T = zeros(m+1,m+1);

% Berechnung der ersten Spalte im Schema
for i = 1:m+1	  % i ist der Zaehler fuer die Stufen
  % Auf jeder Stufe wird die Quadratur ausgefuehrt
  T(i,1) = feval(quadratur, intfunc, a, b, N);
  % von Stufe zu Stufe wird die Anzahl Teilintervalle verdoppelt
  N = 2*N;
end

% eigner source
% T=t(i,j) t(zeilen, spalten)

for j = 2:m+1 %ueber spalten von 2 bis letzter
 for i = 1:m+2-j   % ueber zeilen von 1 bis (m+1)-(j-1) = m+2-j
  T(i,j) = T(i+1,j-1) + 1/(4^j-1) * ( T(i+1,j-1) - T(i,j-1) );
 end
end

result = T(1,m+1);
