function [q,e,nfe] = Quad(f,alpha,beta,nfe)
%  A Gauss-Kronrod (3,7) quadrature over (alpha,beta).

a = [0.2684880898683334, 0.1046562260264672, ...
     0.4013974147759622, 0.4509165386584744];
x = [0.7745966692414834, 0.9604912687080202, 0.4342437493468026];

h = 0.5 * (beta - alpha);
center = alpha + h;
f1 = feval(f,center);
f2 = feval(f,center - h * x(1));
f3 = feval(f,center + h * x(1));
f4 = feval(f,center - h * x(2));
f5 = feval(f,center + h * x(2));
f6 = feval(f,center - h * x(3));
f7 = feval(f,center + h * x(3));
nfe = nfe + 7;

q = h * (5 * (f2 + f3) + 8 * f1) / 9;
qkronrod = h * (a(1) * (f2 + f3) + a(4) * f1 + ...
                a(2) * (f4 + f5) + a(3) * (f6 + f7));
e = qkronrod - q;