function I=trapez(fun,a,b,n)
%I=trapez(fun,a,b,n)
%Computes the integral of the function 
%fun from a to b using the trapese rule
x=linspace(a,b,n);%knots
h=x(2)-x(1);%interval length
fval=feval(fun,x);%function values
w=[1 2*ones(1,n-2) 1];%weights
I=h*w*fval'/2;%trapese formula