
function I=cheb_integral(f,N)

k=(0:N-1)';
x=cos((pi/2)*(2*k+1)/N);
fx=feval(f,x);
I=(pi/N)*sum(fx);