function y = int_rect2(fun,a,b,n)
% numerical integration
% composite rectangle rule
h = (b-a)/n ;
x = a:h:b ;
s = feval(fun,a) ;
for i = 1:n
    s = s + feval(fun,x(i));
end
y = s*h;
end

