function z = polyintcheb(f,a,b)
% POLYINTCHEB uses Chebyshev-spaced 11-point (n=10 degree) polynomial
% interpolation to approximate the integral     /a
%                                               |  f(x) dx
%                                               /b
% This is also known as Clenshaw-Curtis numerical integration.

n = 10;
x = a + ((b-a)/2) * (cos((pi/n)*(n:-1:0)) + 1);
y = feval(f,x);
P = polyfit(x,y,n);        % P(x) approximating f(x)
Q = [P ./ (n+1:-1:1), 0];  % Q'(x) = P(x), and Q has constant term = 0
z = polyval(Q,b) - polyval(Q,a);
