% h2alg
clear all,close all
dx=0.005; L=26;x=0:dx:L; n=length(x);
D2=diag(ones([n-1,1]),-1)-2*diag(ones([n,1]),0)+diag(ones([n-1],1),1);D2=D2/dx^2;
hbar=1;m=0.503985*1822.89;k=1;
alfa=1.0213;D=0.1744;Re=1.3984;
pot=D*(1-exp(-alfa*(x-Re))).^2;
H=(-hbar/(2*m))*D2+diag(pot);[V,E]=eig(H);en=diag(E);
nq=0:n;exa=(0.020*(nq+0.5)-5.7626e-004*(nq+0.5).^2);[en(1:5) exa(1:5)']

for i=1:20
    dE(i)=en(i+1)-en(i);
end

plot(0:17,dE(1:18),'o-')
%axis([0 18 0 0.02])
xlabel('n','FontSize',18)
ylabel('E(n+1)-E(n) / u.a.','FontSize',18)
set(gca,'FontSize',18)
print -deps -r600 SolProb4cap1