function z = polyint(f,a,b)
% POLYINT uses equally-spaced 11-point (n=10 degree) polynomial
% interpolation to approximate the integral     /a
%                                               |  f(x) dx
%                                               /b

n = 10;
dx = (b - a) / n;
x = a:dx:b;
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);
