function demo_normal_domain

a=0; b=2*pi; % Define "normal domain".
psi=@(x) sin(x);
phi=@(x) sin(x)+log(x+3);
f=@(x,y) (x+0.5*y).^10; % integrand
Iex=2.349131532071612e+08; % integral computed with high order rule
n=10; m=11; % Define "n", "m" (cardinality of the rules).

% External routine that computes nodes and weights.
[nodes,weights]=define_normal_rule(n,m,a,b,psi,phi);

% Compute integral.
fnodes=feval(f,nodes(:,1),nodes(:,2));
Inum=weights'*feval(f,nodes(:,1),nodes(:,2));

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

% Plot normal domain (external subroutine)
plot_normal_domain(a,b,psi,phi);
plot(nodes(:,1),nodes(:,2),'go','MarkerEdgeColor','k',...
    'MarkerFaceColor','g','MarkerSize',4);
axis equal; axis tight;
hold off;