function [pt wt]=Rys_pw(c,N,xmin) 
if xmin == 0
     display('xmin = 0; Half range polynomials considered - Enter')
else display('xmin = -1; Full range polynomials considered - Enter')
end
pause
% Quadrature points and weights for Rys polynomials
% xmin = -1 for Full Range J_n polynomials
% xmin = 0 for Half Range R_n polynomials
myfile1=fopen('Rys_alpha_beta.dat','wt')
myfile2=fopen('Rys_pts_wts.dat','wt')
format long e
xmax=1;nint=100;npts=100;xm=[1:1:N];
% pts and wts for the multidomain quadrature
pwmd = mdquad(nint,npts,xmin,xmax); ntot=nint*npts;
p=pwmd(:,1); w=pwmd(:,2); psq=p.*p; 
% weight function w(x)
wtfcn=exp(-c*psq);  
% First two polynomials alpha_0, beta_0 alpha_1, beta_1
b0=ones(ntot,1);
s1=sum(w.*wtfcn); s2=sum(w.*(p.*wtfcn));
h(1)=s1; alfa(1)=s2/s1; beta(1)=0;k=1;
fprintf(myfile1,'Recurrence coefficients, alpha_n and beta_n\n')    
fprintf(myfile1,'%16.12f %16.12f\n',alfa(k),beta(k));
b1=p-alfa(1);
s1=sum(w.*(wtfcn.*(b1.^2))); s2=sum(w.*(p.*(wtfcn.*(b1.^2))));
alfa(2)=s2/s1; h(2)=s1; beta(2)=h(2)/h(1);k=2;
% h(k) are the norms
fprintf(myfile1,'%16.12f %16.12f\n',alfa(k),beta(k));
for k=3:N
    pma=p-alfa(k-1); b2=pma.*b1-beta(k-1)*b0;
    s1=sum(w.*(wtfcn.*(b2.^2))); s2=sum(w.*(p.*(wtfcn.*(b2.^2))));
    alfa(k)=s2/s1; h(k)=s1; beta(k)=h(k)/h(k-1);
fprintf(myfile1,'%16.12f %16.12f\n',alfa(k),beta(k));
    b0=b1; b1=b2;
end
b=beta(2:N); rtb=sqrt(b);
% The Jacobi matrix
J=diag(rtb,-1)+diag(alfa)+diag(rtb,1);
[f,lambda]=eig(J); 
pt=diag(lambda); wt=h(1)*f(1,:).^2; ptwt=[pt,wt']
fprintf(myfile2,'Quadrature points and weights\n')
for i=1:N
    fprintf(myfile2,'%16.12f %1.12e\n',(pt(i)),(wt(i)));
end
% ==============================================================
    function pwmd = mdquad(nint,npts,xmin,xmax)
format long e        
ntot=nint*npts; dx=(xmax-xmin)/ntot;
for i=1:nint
    a=xmin+(i-1)*npts*dx; b=a+npts*dx; pw=fejer2(a,b,npts);
    if i==1
        pwmd=pw;
    else
        pwmd=cat(1,pwmd,pw);
    end
end
% ===============================================================
function pw=fejer2(a,b,N)
format long e
n=N:-1:1; m=1:floor(N/2); th=(2*n-1)*pi./(2*N); p=cos(th');
for k=N:-1:1
  s=sum(cos(2*m*th(k))./(4*(m.^2)-1)); w(k)=2*(1-2*s)/N;
end
r1=(b-a)/2.; r2=(a+b)/2.; ps=r1*p+r2; ws=r1*w; pw=[ps,ws'];