function lagpoly(npoly,alf)
hold off
% Define grid to plot polynomial
xmax=100; ngrid=4001; dx=xmax/(ngrid-1); 
x=[0:dx:xmax]; wtfcn=x.^alf.*exp(-x); srwtfcn=sqrt(wtfcn);
zero2=0*ones(1,ngrid); 
% Define L_0(x_i) and L_1(x_i) and stored as columns in matrix a
a=[]; p1=ones(1,ngrid); a=[x' p1']; p2=(alf+1)*ones(1,ngrid)-x;
a=[a p2'];
% Determine L_n(x_i) by recurrence
for n=3:npoly+2
     nm1=n-1;
    p3=(2*nm1+alf-1-x).*p2/nm1 - (nm1+alf-1)*p1/nm1;
% Store next polynomial in a    
    a=[a p3'];
    p1=p2; p2=p3;
% Redefine p1 and p2 for next iteration in the recurrence relation    
end
% Call to lagptwt for the quadrature points
[pt,wt] = lagptwt(npoly,alf);
npt=length(pt); zero=0*ones(1,npt);
lagpolyn=srwtfcn'.*a(:,npoly+2);
% cmazx is used to scale the y-axis appropriately
cmax=max(lagpolyn);
% Plot sqrt(w(x)L_n(x) versus sqrt(x)
plot(sqrt(x),lagpolyn,'-k','linewidth',1)
hold on
axis([0 10 -0.9*cmax 1.1*cmax])
% Plot quadrature points as circles at the roots of the polynomial
plot(sqrt(pt),zero,'ok','markersize',8,'markerfacecolor','k')
plot(sqrt(x),zero2,'--k','linewidth',2)
xlabel('$\sqrt{x}$','Interpreter','latex','fontsize',20)
ylabel('$\sqrt{x^\alpha e^{-x}}L^{(\alpha)}_n(x)$','Interpreter','latex','fontsize',20)
% Convert N and alpha to strings so as to label the graph
strn1=num2str(npoly);
stra1=num2str(alf);
strn2={'$n = $'};
stra2={'$\alpha =  $'};
% Concatenate the strings
stra=strcat(stra2,stra1);    
strn=strcat(strn2,strn1);
% Add strings to graph
text(7,0.9*cmax,0,strn,'Interpreter','latex','fontsize',20)
text(7,0.7*cmax,0,stra,'Interpreter','latex','fontsize',20)
end
% === Caculates the points and weigths for associated Laguerre polynomials
function [pt,wt] = lagptwt(n,alf)
xn=[0:1:n-1]; a=2*xn+alf+1; b=(xn.*(xn+alf)); rtb=sqrt(b); rtb(1)=[];
t=diag(rtb,-1)+diag(a)+diag(rtb,1);
[f,lambda]=eig(t);
pt=diag(lambda);
for i=1:n
wt(i)=gamma(alf+1)*f(1,i)^2; end
ptwt=[pt,wt'];
end