
function demo_hypercube

ADE=10;
d=3;

% Define integrand
switch d
    case 2
        f=@(x,y) (0.3*x+0.9*y).^10;
        I=5.002201832727280e-01;
    case 3
        f=@(x,y,z) (0.3*x+0.9*y+0.8*z).^10;
        I=4.377443514181815e+01;
end

% Gaussian rule with degree ADE.
n=ceil((ADE+1)/2);
abn=r_jacobi(n,0,0); xw=gauss(n,abn); % Gauss-Legendre
x=xw(:,1); w=xw(:,2);
switch d
    case 2
        [x1,x2] = meshgrid(x); [w1,w2] = meshgrid(w);
        fP=feval(f,x1,x2); w=w1.*w2; 
        Inum=sum(sum(w.*fP));
    case 3
        [x1,x2,x3] = ndgrid(x); [w1,w2,w3] = ndgrid(w);
        fP=feval(f,x1,x2,x3); w=w1.*w2.*w3; 
        Inum=sum(sum(sum(w.*fP)));
end

fprintf('\n \t * I : %1.15e',Inum)
fprintf('\n \t * AE: %1.3e',abs(I-Inum))
fprintf('\n \t * RE: %1.3e \n',abs(I-Inum)/abs(Inum))