
N=100; 
alpha=0.5; beta=1.5;
a=-1; b=1; f=inline('exp(x.^2)');
% a=-1; b=1; f=inline('x.^20');
% a=0; b=1; f=inline('exp(x)');
[x,w]=gauss_jacobi(N,alpha,beta);
xx=a*(1-x)/2+b*(1+x)/2; ww=(b-a)*w/2;
fxx=feval(f,xx);
I=ww'*fxx;
fprintf('\n\t[N]:%3.0f [I]: %1.15e \n\n',N,I)