global c alpha beta
hold on
eps0=.5e-8;
for N=20:5:40
%for N=60:5:80
  for alpha=[-.9:.2:-.1 0 .1:.2:.9 1.5 3 6]
    for beta=[-.9:.2:-.1 0 .1:.2:.9 1.5 3 6]
      if beta>=alpha
        for i=1:9
          c=i/10;
          ab=r_subjacobi(N,eps0,c,alpha,beta); 
%          ab=r_subjacobi0(N,eps0,c,alpha,beta);
          xw=gauss(N,ab);
          x=xw(:,1); w=xw(:,2);
          y=N*w./(pi*(((1-c*x).^alpha).*((1+c*x).^beta)));
%           y=N*(2/(1+c))^(alpha+beta)*w./(pi ...
%       *((3-c)/(1+c)-x).^alpha.*(1+x).^beta);
          plot(x,y,'.')
          plot([-1.1 1.1],[0 0])
          plot([0 0],[0 1],'--')
          axis('equal')
          axis([-1.1 1.1 -.1 1.2])
        end
      end
    end
  end
end
