function [I,x,w]=integration_gauss_jacobi(N,alpha,beta,f)
% INPUT:
% N: number of nodes.
% alpha,beta: w(x)=((1-x)^alpha)*((1+x)^beta)
% f: function such that the integrand is f(x) w(x) in (-1,1).
% OUTPUT:
% I: approximation of I(f*w,a,b)
% x,w: nodes and weights of the rule.
ab=r_jacobi(N,alpha,beta); % recursive terms
xw=gauss(N,ab); % matrix of nodes and weights
x=xw(:,1); w=xw(:,2);% nodes and weights
fx=feval(f,x); % function evaluation
I=w'*fx; % value of the integral
