% to carry out a variational calculation of the energy levels of the
% quartic oscillator

% first two hermite polynomials; 
% note that in atomic units alph=mu*omega = sqrt(k*mu)

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:6;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^2 dx) = int(psi*psi*u^2 du)/(alph^3/2)
disp('calculating matrix of x^2 ...')
x2mat=zeros(6,6);
% 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
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
% these are matrix elements of u^2.  we need to convert to matrix elements
% of (R-Re)^2.

x2mat=x2mat/alph;
disp('done')
% now determine matrix elements of x
disp('calculating matrix of x^1 ...')
x1mat=zeros(6,6);
syms x1mat
for i=1:6
   for j=i:1:6
      x1mat(i,j)=int(psi(i).'*psi(j)*u,u,-inf,inf)/(alph^(3/2));
      x1mat(j,i)=x1mat(i,j);
   end
end
x1mat=x1mat/sqrt(alph);
% now determine matrix elements of x^3
disp('done')
disp('calculating matrix elements of x^3 ...')
x3mat=zeros(6,6);
syms x3mat
for i=1:6
   for j=i:1:6
      x3mat(i,j)=int(psi(i).'*psi(j)*u*u*u,u,-inf,inf)/(alph^(3/2));
      x3mat(j,i)=x3mat(i,j);
   end
end
disp('done')
% now determine matrix elements of -d^2/dx^2
% we know that d/dx= du/dx * d/du = sqrt(alph) d/du
% so that d^2/dx^2= alph * d^2/du^2

% also 
% int(psi d^2/dx^2 psi dx) = (alph/sqrt(alph)) int(psi d^2/du^x psi du)
disp('calculating matrix elements of -d^2/dx^2')
for i=1:6
   for j=i:2:6
      k2mat(i,j)=int(psi(i).'*diff(psi(j),u,2),u,-inf,inf)*sqrt(alph);
      k2mat(j,i)=k2mat(i,j);
   end
end

% take into account that the kinetic energy operator is -(d^2/dx^2)/(2*mu)

syms mu k2

k2mat=-k2mat/(2*mu);
disp('done')

% to check that the matrix of the harmonic oscillator hamiltonian is
% diagonal in the basis of the harmonic oscillator functions

h2mat=k2mat+x2mat*k2/2; h2mat=simplify(h2mat);
h2mat=subs(h2mat,alph,sqrt(k2*mu));h2mat=simplify(h2mat);

% this is diagonal with values of (n+1/2)sqrt(k2/mu) = (n+1/2)omega
% which is what we would expect
disp('calculating matrix of x^4 ...')      
% 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
disp('done')  
syms  k4

% matrix of quartic oscillator hamiltonian is
h4mat=k2mat+x4mat*k4/2; h4mat=simplify(h4mat);

% now, apply perturbation theory, where V=0.5k4x^4 - 0.5k x^2

% we just need the first row of the matrix, call this v2.  these are the
% matrix elements of H' between the harmonic oscillator functions

v2=h4mat(1,:)-h2mat(1,:);
% determine numerical values of this
vv2=subs(v2,[k2 k4 mu],[0.37 5.56 918]);
% and substitute for alpha=sqrt(mu*k2)
vv2=subs(vv2,alph,sqrt(0.37*918));
% E(0) is the harmonic oscillator energy
E0=0.5*sqrt(.37/918);

% E(1) is just the first term in vv2
E1=vv2(1);

%E(2) is the square of the other terms divided by their energy denominators
edenom=-[1:4]*sqrt(.37/918); % note negative sign: E(0)-E(n)
E2=sum(vv2(2:5).^2./edenom);

% disp('energy of quartic oscillator through 2nd order
disp('     E0           E1           E2          E0+E1+E2')
disp([E0,E1,E2,E0+E1+E2])


%now for variational calculation


% determine numerical matrix of quartic oscillator hamiltonian when
%mu=918;k4=5.56;k2=0.37  

% and alph=sqrt(mu*k2)

hmat=subs(h4mat,[k4 mu ],[5.56 918 ]);

hmat=subs(hmat,alph,sqrt(0.37*918));

[evec eval]=eig(hmat)
% the matrix evec contains the eigenvectors (by column)
% the diagonal matrix eval is the matrix of eigenvalues

% note that evec'*h4mat*evec = eval

% to plot out the wavefunction and compare it with the H.O. wavefunction
syms x
v=0;
%this is the harmonic oscillator wavefunction
psiho=subs(psi(v+1),u,sqrt(alph)*x);
psiho=subs(psiho,alph,sqrt(0.37*918));
% for the quartic oscillator

psiqo=subs(psi*evec(:,v+1),u,sqrt(alph)*x);
psiqo=subs(psiqo,alph,sqrt(0.37*918));
% this plots the harmonic (dashed curve) and quartic oscillator (solid curve)
% wavefunctions
ezplot(psiqo);hold;h=ezplot(psiho);set(h,'linestyle','--');hold


% now allow alpha to vary.  note that we don't try to do this symbolically,
% since we can't obtain a symbolic representation of the roots of a 6 x 6
% matrix.  We do it numerically, and then fit the results in the region of
% the minimum, and then differentiate this polynomial fit

alpha=[18:2:60];
for i=1:length(alpha)
   eval=eig(subs(h4mat,[k4 mu alph],[5.56 918 alpha(i)]));
   ee(i)=eval(1);
end

figure(2);
plot(alpha,ee);
set(gca,'fontsize',12);xlabel('alpha');ylabel('Energy');
title('5 hermite function variational approx to quartic oscillator')

% determine the minimum on the grid of alpha that we have used
[e mn]=min(ee);
% fit E(alpha) to a quadratic
cc=polyfit(alpha(mn-1:mn+1),ee(mn-1:mn+1),2);
% find the roots of dE(alpha)/dalpha
alph_min=roots(polyder(cc));
% evaluate the value of E(alph_min) using the previously determined
% polynomial fit.
emin=polyval(cc,alph_min);

% check by using alph_min in the full calculation
hmat=subs(h4mat,[k4 mu alph],[5.56 918 alph_min]);
eval=eig(hmat);

disp(['lowest energy = ' num2str(min(eval)) ', obtained for alpha= ' num2str(alph_min)])

% determine and plot wavefunction
[evec eval]=eig(hmat);
psiqqo=subs(psi*evec(:,v+1),u,sqrt(alph)*x);
psiqqo=subs(psiqqo,alph,alph_min);
figure(1)
ezplot(psiqo);hold;h=ezplot(-psiqqo);set(h,'color','m');h=ezplot(psiho);set(h,'linestyle','--');hold
set(gca,'fontsize',12);title(['blue: \alpha=' num2str(sqrt(0.37*918)) '; magenta: \alpha=' num2str(alph_min) '; dashed: harmonic oscillator'])


root=fzero('0.37*x^2-5.56*x^4',-.2);
hold;
plot(root*[1 1],[0 1.4],'-k');hold


