
%calculates, using matlab's symbolic package, the matrix elements of x^2
%and x^4 in a harmonic oscillator basis
% 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:9;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)

clear x2mat
% 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:6
   for j=i:2:6
      x2mat(i,j)=int(psi(i).'*psi(j)*u*u,u,-inf,inf)/(alph^(3/2));
      x2mat(j,i)=x2mat(i,j);
   end
end
      
% calculate matrix of x^4
for i=1:6
   for j=i:2:6
      x4mat(i,j)=int(psi(i).'*psi(j)*u^4,u,-inf,inf)/(alph^(5/2));
      x4mat(j,i)=x4mat(i,j);
   end
end
   
