function ShowGK()
% Illustrates the 15-7 Gauss-Kronrod Pair
a = 2;
b = 5;
[QK,QG] = GaussKronrod(@f,a,b)
I = quad(@f,a,b,.000000000001);
clc
fprintf('Integrate f(x) = (1./(1+3*x.^4)).*cos(x).*x from 2 to 5\n\n')
fprintf('QK = %20.15f\nQG = %20.15f\nQ = %20.15f\n\n',QK,QG,I)
err = 200*abs(QG-QK)^1.5;
fprintf('Heuristic error estimate = %10.3e\n',err)
fprintf('Actual error using quad = %10.3e\n',abs(QK-I))

function y = f(x)
y = (1./(1+3*x.^4)).*cos(x).*x;

function [QK,QG] = GaussKronrod(f,a,b)
% The 15 nodes...
xGK = [...
-0.991455371120813;...
-0.949107912342759;...
-0.864864423359769;...
-0.741531185599394;...
-0.586087235467691;...
-0.405845151377397;...
-0.207784955007898;...
0.000000000000000;...
0.207784955007898;...
0.405845151377397;...
0.586087235467691;...
0.741531185599394;...
0.864864423359769;...
0.949107912342759;...
0.991455371120813];

% The Kronrod weights...
wK = [...
0.022935322010529;...
0.063092092629979;...
0.104790010322250;...
0.140653259715525;...
0.169004726639267;...
0.190350578064785;...
0.204432940075298;...
0.209482141084728;...
0.204432940075298;...
0.190350578064785;...
0.169004726639267;...
0.140653259715525;...
0.104790010322250;...
0.063092092629979;...
0.022935322010529];

% The Gauss weights associated with x(2:2:14)...
wG = [...
0.129484966168870;...
0.279705391489277;...
0.381830050505119;...
0.417959183673469;...
0.381830050505119;...
0.279705391489277;...
0.129484966168870];

x = ((a+b)/2)+((b-a)/2)*xGK
fvals = f(x);
QK = (wK'*fvals)*(b-a)/2;
QG = (wG'*fvals(2:2:14))*(b-a)/2;