function bimode_polyn_with_pts(npoly)
hold off
eps=0.005;
xmax=2; ngrid=1501; dx=2*xmax/(ngrid-1); x=[-xmax:dx:xmax]; 
xsq=x.*x; x4=xsq.*xsq; wtfcn=exp(-(x4/4-xsq/2)/eps);srwtfcn=sqrt(wtfcn);
%grid to plot polynomial and root of the weight function
load bimodal_eps005.dat; n=100;           %delete the bottom m+1 to n rows
b=bimodal_eps005(:,2);
g=bimodal_eps005(:,3);
rtb=sqrt(b);  
rtbx=rtb;
rtbx(1)=[];
rtbx(npoly:n-1)=[];
%delete last element of the off-diagonal elements
t=diag(rtbx,-1)+diag(rtbx,1);
[f,lambda]=eig(t);
pt=diag(lambda);
%pause
for i=1:npoly
wt(i)=g(1)*f(1,i)^2; end
ptwt=[pt,wt'];
zero2=0*ones(1,ngrid); %To draw dashed line at y = 0 
a=[];
gam=sqrt(pi/4); srgam=sqrt(gam); rtpi=sqrt(pi);
%Create the first polynomial as a vector with unit components:
p1=ones(1,ngrid)/sqrt(g(1));
%Create the matrix A with the first column x and the second p(1):
a=[a p1'];
%Create the second polynomial as a vector:
norm2=ones(1,ngrid)/sqrt(g(2));
p2=x/sqrt(g(2));
%Add this vector to the 3rd column of A:
a=[a p2'];
for n=2:npoly
%Use the recursion relation to get the next polynomial:
    p3=((x.*p2)- p1*rtb(n))/rtb(n+1);
    %Add this polynomial to the next column of A:    
    a=[a p3'];
%Store the previous two polynomials so as to use in the recursion
%the next time around in the loop
    p1=p2; p2=p3;
end
npt=length(pt); zero=0*ones(1,npt);
%maxpolyn=srwtfcn'.*a(:,npoly+2);
maxpolyn=srwtfcn'.*p3';
plot(x,maxpolyn,'-k','linewidth',1.6)
hold on
plot(pt,zero,'ok','markersize',7,'markerfacecolor','k')
hold on
plot(x,zero2,'--k','linewidth',1.2)
axis([-xmax xmax -2 2])
set(gca,'FontSize',36)
set(gca,'Ytick',[-2:.5:2],'linewidth',1.6)
set(gca,'Xtick',[-xmax:1:xmax],'linewidth',1.6)
xlabel('$x$','Interpreter','LaTex','FontSize',36)
ylabel('$\sqrt{w(x)}B_n(x)$','Interpreter','LaTex','FontSize',36)
strn1=num2str(npoly);
strn2={'$n = $'};    
strn3=strcat(strn2,strn1);
text(-.5,-1.7,strn3,'Interpreter','Latex','Fontsize',36)
strn4 = {'$\epsilon = 0.005$'};
text(-.5,1.7,strn4,'Interpreter','Latex','Fontsize',36)
set(gcf, 'Units','centimeters','Papersize',[36,36])
set(gcf, 'Units','centimeters','Position',[3 3 24 20])