% first two hermite polynomials
syms h u alph
h(1)=1;
h(2)=2*u;
% recursion relation:  H_{n+1}=2[u H_n - nH_{n-1}
for i=3:5;h(i)=2*(u*h(i-1)-(i-2)*h(i-2));end


h=simplify(h);


s2=int(h.*h*exp(-u^2),-inf,inf);
%normalization constant
n=1./sqrt(s2);n=simplify(n);


psi=n.*h*exp(-0.5*u^2);

% now, we want int(psi^2 dx) to equal one, not int(psi^2 du)
% where u=sqrt(alph)*x

% since du=sqrt(alph)dx,  we need to multiply the wavefunctions
% by (alph)^(1/4) to be normalized with respect to integration over x

psi=psi*alph^(1/4);

% thus int(psi*psi dx)= int(psi*psi du)/sqrt(alph) = 1

% evaluate matrix of x^2
% x^2 dx = alph^(3/2  u^2 du, 
% so int(psi*psi*x*x dx) = int(psi*psi*u*u du)/(alph^3/2)
 
x2mat=zeros(5,5);
% this takes into account that by symmetry only even-even or odd-odd
% elements are nonzero. also we explicity calculate only the upper triangle
% of the matrix
syms x2mat x4mat
for i=1:5
   for j=i:2:5
      x2mat(i,j)=int(psi(i).'*psi(j)*u*u,u,-inf,inf)/(alph^(3/2));
      x2mat(j,i)=x2mat(i,j);
   end
end
      
% and similarly for matrix of x^4
for i=1:5
   for j=i:2:5
      x4mat(i,j)=int(psi(i).'*psi(j)*u^4,u,-inf,inf)/(alph^(5/2));
      x4mat(j,i)=x4mat(i,j);
   end
end
   

% to plot wavefunctions

wf=subs(psi(3),alph,1)

ezplot(wf)

% to calculate the term <v=0|H'|v=n>^2/(E_0-E_n), which appears in the
% expression for the perturbation correction to 2nd order

mu=918;k2=0.37;alpha=sqrt(k2*mu);

k4=5.56;

% let vvec be the vector <v=0|H'|v=n>

vvec=0.5*(x4mat(:,1)*k4-x2mat(:,1)*k2)

% substitute for alph the particular value alpha
vvec=subs(vvec,alph,alpha)

% turn this into a row vector
vvec=vvec';

% square and divide by energy differences

% energies of harmonic oscillator states are En=(n+1/2)*hbar*omega= (in
% atomic units) (n+1/2)*sqrt(k/mu) 
e=([0:4]+1/2)*sqrt(k2/mu);

% energy denominators (En-E0)
e_denom=e(1)-e(2:5);

vvec(2:5).^2./e_denom

% so first order correction

E1=vvec(1);

% 2nd order correction

E2=sum(vvec(2:5).^2./e_denom);

% convert to cm-1 units

[E1 E2]*2194784.6;


