function [x,w]=radau(al,be,N,x_zero,x_uno,x,p,d1p)

%  W = quadratura pesos
%     1 = quadratura de radau ( incluje x=1)
%     2 = quadratura de radau (incluje x=0)


Npts=length(x);
% RADAU  incluje (x=1);    quadfraure ==1

if ((x_zero==0)&(x_uno~=0))
    for i=1:Npts
        w(i,1)=(1/x(i))/(d1p(i)^2);
    end
    w(Npts,1)=w(Npts)/(al+1);
end  % end de if

% RADAU  incluje (x=0);    quadraure ==2

if ((x_zero~=0)&(x_uno~=0))
    for i=1:Npts
        w(i,1)=(1/(1-x(i)))/(d1p(i)^2);
    end
    w(1,1)=w(1)/(be+1);
end  % end de if


% RADAU  incluje (x=0) & (x=1);    quadraure ==3

if ((x_zero~=0)&(x_uno~=0))
    for i=1:Npts
        w(i,1)=d1p(i)^(-2);
    end
    w(1,1)=w(1)/(be+1);
    w(Npts,1)=w(Npts)/(a1+1);
end  % end de if
    
% Equ 8.10 b com gamma(N,N) dada na equa  8.85 a
w=(gamma(al+1)*gamma(be+1)/gamma(al+be+2))*w/sum(w);
end