function demo_kronrod(f,alpha,beta,N)

% Integration of a function 'f' in [-1,1] and Gauss-Kronrod rules.
if nargin < 1, f=@(x) exp(x).*sqrt(1-x.^2); end
if nargin < 3, alpha=0; beta=0; end % Jacobi weight
if nargin < 4, N=10; end

dim=ceil(3*N/2)+1; % length Kronrod rule
ab=r_jacobi(dim,alpha,beta); % Gauss-Legendre recurrence (as required by kronrod)
xwg=gauss(N,ab); xwk=kronrod(N,ab);

g=@(x) f(x).*(1-x).^alpha.*(1+x).^beta;
I=integral(g,-1,1,'AbsTol',10^(-15),'RelTol',10^(-15));

fk=feval(f,xwk(:,1)); % this evaluation includes that at Gauss-Legendre nodes!

Sg=(xwg(:,2))'*fk(2:2:end);
Sk=(xwk(:,2))'*fk;

fprintf('\n * Gauss-Legendre rule \n \n'); xwg
fprintf('\n * Gauss-Legendre-Kronrod rule \n \n'); xwk
fprintf('\n \t I gauss : %1.15e',Sg);
fprintf('\n \t I kronr : %1.15e',Sk);
fprintf('\n \t I exact : %1.15e',I);
fprintf('\n \n');
fprintf('\n \t E gauss : %1.3e',abs(I-Sg));
fprintf('\n \t E estim : %1.3e',abs(Sk-Sg));
fprintf('\n \n');
