function Harm_Osc_Herm2(nmax)
% Exact calculation of the eigenvalues and eigenfunctions of the 
% quantum harmonic osciallator with Hermite polynomials; Eq. (6.178)
% With nmax = 12,  this code reproduces Fig. 6.15
% This code could be vectorized and shortened.
hold off
format long e
npoly=nmax; rtpi=sqrt(pi);
[ptwt] = hermptwt(nmax);
pt = ptwt(:,1); wt=ptwt(:,2); omptsq=diag(1-pt.*pt);
% Store H_n(x_i) and H'_n(x_i), n=0 and n=1 in matrices poly and ployp
xn=[0:1:nmax]; b=xn/2; rtb=sqrt(b); v=ones(1,nmax); poly=[]; polyp=[]; 
poly=[poly (v./sqrt(sqrt(pi)))']; poly=[poly (pt*sqrt(2/sqrt(pi)))];
w = zeros(1,nmax); polyp=[polyp w']; polyp=[polyp (v.*sqrt(2/sqrt(pi)))'];
% Calculate H_n(x_i) and H'_n(x_i) by recurrence
% Code could be rewritten to use function [hn] below
for n=2:nmax-1
    xp=pt.*poly(:,n)/rtb(n+1)-poly(:,n-1)*rtb(n)/rtb(n+1);
        eigf1=xp.*exp(-pt.*pt/2); eigf3=xp;
    poly=[poly xp];
    yp=pt.*polyp(:,n)/rtb(n+1)-polyp(:,n-1)*rtb(n)/rtb(n+1)+poly(:,n)/rtb(n+1);
    polyp=[polyp yp];
end
% Calculate the symmetrized derivative operator as given by Eq. (3.138)
for i=1:nmax
    xppi=polyp(i,:);
    for j=1:nmax
        ypj=poly(j,:); s=sqrt(wt(i)*wt(j))*sum(ypj.*xppi);
d(i,j)=s;
    end
end
% mat2 is the Hamiltonian, Eq. (6.178).
mat2=d'*d;
[eigfcns,eigen]=eig(mat2);
eig1=diag(eigen)';
eig2=sort(eig1,'ascend');
display(' n     Approx      Approx/Exact')
for i=2:nmax
   fprintf('%2i %12.8f %12.8f\n' ,i-1,eig2(i),eig2(i)/(2*(i-1)));
end
eigf2=exp(-pt.^2/2).*eigfcns(:,nmax)./sqrt(wt);
xnorm1=sum(wt.*eigf1.^2); xnorm2=sum(wt.*eigf2.^2);
factor=sqrt(xnorm1/xnorm2);
[hn]= hermite(npoly+1);
xmax=10; ngrid=401; dx=2*xmax/(ngrid-1); x=[-xmax:dx:xmax];
% Plot the Hermite polynomial
plot(x,hn(:,npoly),'-k','linewidth',1.)
hold on
% Plot the corresponding eigenfunctions
plot(pt,eigf2,'sk','linewidth',1.6,'markersize',10,'markerfacecolor','k')
% Format the graph
axis([-8 8 -.6 .6])
set(gca,'FontSize',36)
set(gca,'Ytick',[-.6:.2:.6],'linewidth',1.6)
set(gca,'Xtick',[-8:2:8],'linewidth',1.6)
xlabel('$x$','Interpreter','latex','fontsize',36)
str2=num2str(nmax);
ylabel('${\rm H}_{{\rm n}}(x)e^{-x^2/2}$','Interpreter','latex','fontsize',36)
str1={'n = '};
str3=strcat(str1,str2);
text(2,-.5,-1,str3,'Interpreter','latex','fontsize',36)
set(gcf, 'Units','centimeters','Papersize',[36,36])
set(gcf, 'Units','centimeters','Position',[3 3 24 20]) 
% ===================== HERMITE QUAD PTS and WTS ================
function [ptwt] = hermptwt(n)
xn=[0:1:n-1]; b=xn/2; rtb=sqrt(b); rtb(1)=[]; t=diag(rtb,-1)+diag(rtb,1);
[f,lambda]=eig(t); pt=diag(lambda); wt=sqrt(pi)*f(1,:).^2;
ptwt=[pt,wt'];

%================== ORTHONORMAL HERMITE POLYNOMIALS ======================
function [hn]= hermite(npoly)
format long e
xn=[0:1:npoly]; b=xn/2; rtb=sqrt(b);
xmax=10; ngrid=401; dx=2*xmax/(ngrid-1); 
x=[-xmax:dx:xmax]; wtfcn=exp(-x.*x); srwtfcn=sqrt(wtfcn);
v=ones(1,ngrid);
hn=[]; hn=[hn (v.*srwtfcn/sqrt(sqrt(pi)))'];
hn=[hn (x.*srwtfcn*sqrt(2/sqrt(pi)))'];
for n=2:npoly
%Recurrence relation 
pnext=x'.*hn(:,n)/rtb(n+1)-hn(:,n-1)*rtb(n)/rtb(n+1);
    hn=[hn pnext];
end