% -------------------------------------------
%
% Romberg für das Testbeispiel
% HeSe 2012/13
%
% ----------------------------------------

clear all

f = inline('sin(x)')      % a = 0 und b = pi/2
%f = inline('1/x')                   % a = 1 und b = 2

a   = input('untere Intervallgrenze a = ');
b   = input('obere  Intervallgrenze b = ');
eps = input('Toleranz eps = ');

Nmax = input('maximale Anzahl Halbierungen Nmax = ');

k = 2;

h = b - a;              % h0
T0 = h/2*(f(a) + f(b));
M0 = h*f(a + h/2);
T1 = 1/2*(T0 + M0);
h = h/2;                % h1
n = 2;
R(1,1) = T0;
R(2,1) = T1;
R(2,2) = 1/3*(4*R(2,1) - R(1,1));
quot = [];
while abs(R(k,k) - R(k-1,k-1)) > eps & k < Nmax
    k = k + 1;
    h = h/2;            % h2
    n = 2*n;
    mitten = 0.0;
    for i = 1:2:n
         mitten = mitten + f(a+i*h);
    end
    mitten = 2*h*mitten;
    R(k,1) = 1/2*(R(k-1,1) + mitten);
    
    % Romberg Schema
    p = 1;
    for j = 2:k
        p = 4*p;
        R(k,j) = (p*R(k,j-1) - R(k-1,j-1))/(p-1);
        
    end
    
   
end

