function Int=IRQ(w,x,f)
%-----------------------------------------------------------------
% This function computes the integral of f with Indirect Rational
% Quadrature (IRQ method).
% 
% INPUTS:
% x - vector of nodes
% w - vector of F-H's weights
% f - integrand function
% 
% OUTPUTS:
% Int - integral of f
% 
%------------------------------------------------------------------
% Authors: 
%      Stefano De Marchi and Cinzia Bandiziol
%      Department of Pure and Applied Math.
%      University of Padova
%   
%           
% Version: May 2015		
%------------------------------------------------------------------	

[a, b]=size(x);
	if (a==1 && b>1)
		x=x';
    end
    
n=size(x,1);
% Compute the constant term
constant_term_old=feval(f,x);
constant_term=constant_term_old(2:end);
T_old=repmat(x,1,n)-repmat(x',n,1);
	
T_new=ones(n,n)./T_old;
	
[a, b]=size(w);
	if (a==1 && b>1)
		w=w';
	end
c=find(w==0);
	if (size(c,1)>=1)
		printf('A weight is null');
		return
	end			
W1=repmat(w',n,1);% nxn matrix
W2=repmat(w,1,n);% nxn matrix
W=W1./W2;
D_old=T_new.*W;

a=find(D_old==Inf);
D_old(a)=0;
D_vect=(-1)*sum(D_old,2);

% Compute a matrix with D_vect as diagonal
D=spdiags(D_vect,0,D_old);
	
D_new=D(2:end,2:end);

% Compute the approximation of antiderivative of f in x
u=D_new\constant_term;
Int=u(n-1);