% Berechnung der
% Trapezwerte


%  Def der Funktion im m-file f.m

format long

% f = inline('exp(-x^2)')
f = inline('sin(x)')

tt = clock; %flops(0);
 
a = 0;			% untere Grenze
b = 2*atan(1);    % obere Grenze
%b = 0.5

eps = 1e-5;		% Toleranz
kmax = 15;		% max Anzahl Halbierungen

h= b-a;

T = [];
M=0;
t = h*(f(a) + f(b))/2; 	% 0-ter Trapezwert

T = [T,t];

k = 1;
n = 1;
  
fprintf('\n');
while  k<=kmax & abs(t-M)>eps*(t+1)	% Mischung aus
					% relativem und
					% absolutem Fehler
  k = k+1;  
  h = h/2;	% neue Pommes-Frites-Dicke
  n = 2*n;

  % Berechnung der neuen Mitten
  M = 0;
  for j= 1:n/2				% n/2=Anz neue Mitten
    M = M + f(a+(2*j-1)*h);
  end
  M = 2*h*M;				% da bereits mit der
					% neuen pommes-frites
					% Dicke gerechnet wird
  
  t = (T(k-1) + M )/2;			% neuer Trapezwert
  fprintf('%5.0f %5.0f %15.8f %20.15f \n',k-1, n, h ,t);
  T = [T,t];
end
  
%fprintf('ntime = %4.2f secs   flops = %6.0f\n', etime(clock,tt),flops);
fprintf('ntime = %4.2f secs   \n', etime(clock,tt));