function f = rechteck( fun, a, b )
    N = 1000;            % 1000 Teilintervalle
    d = (b-a)/N;         % Breite der Teilintervalle
    x = a:d:b;           % Vektor mit den Stützpunkten
    f = 0;               % Startwert der Summation
    % Summation über alle Intervalle
    for( n=1:N ) 
       f = f + fun(x(n)) * d;  % Funktionswert * Breite
    end
end