% to determine wavefunction and energies for
% quartic oscillator using discrete (grid) representation
% of the wavefunction
% define mass
% H2 reduced mass (0.5 amu)
  mu=918.58;
% equipositioned grid extending delr bohr on each side of
% the minimum
%  delr=0.6;  % assume delr is predefined
% equilibrium internuclear distance for H2
re=1.4011;
% delr defines the range of integration
% re-delr< r < re+delr
delr=0.8;
rmin=re-delr;
rmax=re+delr;
% hh is the approximate spacing between points
hh=0.01;  
% determine approximate number of points
% the 'floor' function converts M to an integer
M=floor(2*delr/hh)+1;

% readjust h slightly so that the range is preserved
h=(rmax-rmin)/(M-1);
r=[rmin:h:rmax];



% potential for quartic oscillator
% force constant is 
k4=2.8417;
v=0.5*k4*(r-re).^4;
% construct diagonal matrix with v(r) along diagonal
pemat=diag(v);
% determine kinetic energy matrix using 3 point differentiation
% formula  d2f/dx^2 = [f(x-1) -2 f(x) + f(x+1)]/h^2
kemat=zeros(M,M);
kemat=diag(ones(M-1,1),1);
kemat=kemat+kemat';
kemat=kemat-2*diag(ones(M,1),0);
kemat=kemat/(h^2);
kemat=-kemat/(2*mu);

% % here are the equivalent statements for the 5 point approximation
% kemat=zeros(M,M);
% kemat=16*diag(ones(M-1,1),1);  
% kemat=kemat-diag(ones(M-2,1),2);
% kemat=kemat+kemat';
% kemat=kemat-30*diag(ones(M,1),0);
% kemat=kemat/(12*h^2);
% kemat=-kemat/(2*mu);
% % 

% total hamiltonian is sum of kemat and pemat
hmat=kemat+pemat;
% determine eigenvalues and eigenvectors directly of hmat
% using matlab call "eig"

[evec eval]=eig(hmat);

% normalize eigenvectors (sum of squares equals 1 already)
n=1/sqrt(h);
evec=evec*n;
% sort eigenvalues
[energy nn]=sort(diag(eval));
% first two energies
format shorte
[h M hh energy(1:2)']
format short
% plot wavefunction if requested
  plot(r,evec(:,1:2),'-o');

   

