function Harm_Osc_Herm1
% Scan four values of nmax 6, 7, 8 and 9 as in Table 6.6
for nmax=6:1:9
fprintf(' nmax = %2i\n',nmax)
% Hermite quadrature points and weights
[ptwt]=hermptwt(nmax);
pt=ptwt(:,1);
xn=[1:1:nmax];
% Hamiltonian for harmonic oscillator as in Eq. (6.175)
for i=1:nmax
    exact(i)=(i-0.5);
    for j=1:i
        if i==j 
            h(i,j)=(4*nmax-1-2*pt(i)^2)/12+pt(i)^2/2;
        else
            h(i,j)=(-1)^(i-j)*(1/(pt(i)-pt(j))^2-0.25);
        end
        h(j,i)=h(i,j);
    end
end
% Eigenvalues
[f,en]=eig(h); lambda=diag(en);
display('n      Exact     Approx')
for i=1:nmax
fprintf('%i %10.5f %10.5f\n',i,exact(i),lambda(i));
end
end
% ================ 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'];