function Harm_Osc_Fourier
clc
% Choose one of 4 dx values
dxvec=[.5 .25 .125 .0625];
% String synmbols to label curves
str1='-ok'; str2='-sk'; str3='-^k'; str4='-ok';
symbol=[str1 str2 str3 str4]; xmax=7; fac=pi^2/3;
kmax=0;
% Scan through the four dx values
for kk=1:4
    kmin=kmax+1; kmax=kmin+2;
    dx=dxvec(kk); x=[-xmax:dx:xmax]; nmax=length(x); xn=[1:1:nmax];
% Hamiltonian for harmonic oscillator with Fourier basis functions
for i=1:nmax
    exact(i)=(i-0.5);
    for j=1:nmax
        if i==j 
            h(i,j)=fac/(2*dx^2)+x(i)^2/2;
        else
            h(i,j)=(-1)^(i-j)/(dx^2*(i-j)^2);
        end
    end
end
% eigenvalues
[f,en]=eig(h);
lambda=diag(en);
for i=1:nmax
acc(i)=log10(abs(1-lambda(i)/exact(i)));
%fprintf('%i %16.8e %16.8e %16.8e\n',i,exact(i),lambda(i),acc(i));
end
% Plot the accuracy versus N to the 4 dx values
nm=20;
if mod(kk,2) == 0
plot(xn(1:nm),acc(1:nm),symbol(kmin:kmax),'markersize',8)
hold on
else
plot(xn(1:nm),acc(1:nm),symbol(kmin:kmax),'markersize',8,'markerfacecolor','k')
end
hold on
xlabel('${\rm n}$','Interpreter','latex','fontsize',32)
ylabel('$\log_{10}[{\rm Relative}\;\; {\rm Error}]$','Interpreter','latex','fontsize',32)
axis([0 20 -16 0])
set(gca,'FontSize',32)
set(gca,'Ytick',[-16:4:0],'linewidth',1.6)
set(gca,'Xtick',[0:5:20],'linewidth',1.6)
str={'(A)'};
text(16,-1,str,'Interpreter','latex','fontsize',32)
leg=legend('$N=28$','$N=56$','$N=112$','$N=224$','Location','NorthWest');
set(leg,'Interpreter','latex','fontsize',26);
str={'$\alpha = 4; s=1/\sqrt{\alpha} = 0.5$'};
text(.18,6,str,'Interpreter','latex','fontsize',32)
set(gcf, 'Units','centimeters','Papersize',[36,36])
set(gcf, 'Units','centimeters','Position',[3 3 24 20]) 
end
% ===================== HERMITE QUAD PTS and WTS ================
function [ptwt] = hermptwt(n)
xn=[0:1:n-1];
b=xn/2; rtb=sqrt(b);
rtb(1)=[];
rtb;
t=diag(rtb,-1)+diag(rtb,1);
%disp(t);
[f,lambda]=eig(t);
pt=diag(lambda);
for i=1:n
wt(i)=sqrt(pi)*f(1,i)^2; end
ptwt=[pt,wt'];
s=sum(wt);