function y = int_trap(fun,a,b,n)
% numerical integration
% composite trapezium rule
h = (b-a)/n ;
x = a:h:b ;
s = 0 ;
for i = 1:n
    s = s + feval(fun,x(i));
end
fa = feval(fun,a);
fb = feval(fun,b);
y = ((fa+fb)/2+s)*h;
end

