function [Q,N,X,quad,err] = gauss(f,a,b,A,W,toler)
% [Q,N,X,quad,err] = gauss(f,a,b,A,W,toler)
% Gaussian quadrature.
% f is the name of the function, input.
% a is the left endpoint, input.
% b is the right endpoint, input.
% A is the table of abscissas, input.
% W is the table of weights, input.
% toler is the tolerance, input.
% Q is the table of values, output.
% N is the number of abscissas used, output.
% X is the abscissas that were used, output.
% quad is the quadrature value, output.
% err is the error estimate, output.
mid = (a+b)/2;
wide = (b-a)/2;
err = 1;
j = 1;
x = mid + A(1,1);
N(1) = 1;
Q(1)= W(1,1)*feval(f,x)*wide;
while (((err>toler)|(j<4))&(j<17))
j = j+1;
sum1 = 0;
i = j;
if (j>10),
i = 12 + 4*(j-11);
end
if (j>14),
i = 24 + 8*(j-14);
end
for k = 1:i,
X(k) = mid + A(j,k)*wide;
sum1 = sum1 + W(j,k)*feval(f,X(k));
end
N(j) = i;
Q(j) = sum1*wide;
quad = Q(j);
err = abs(Q(j)-Q(j-1));
end