%Disclaimer: zero efforts have been made to optimize anything in the notebook below 
%cutoff parameter
nn = 10;
%dimensionless interaction strength 
mu = 0.1;
%Hamiltonian Matrix
ham = zeros(nn); 
for n = 0:nn
    for m = 0:nn
        ham(n+1,m+1)=( sqrt(n*(n-1)*(n-2)*(n-3))*kronecker(m,n-4) + (4*n-2)*sqrt(n*(n-1))*kronecker(m,n-2) + (6*n^2 + 6*n + 3)*kronecker(m,n) + sqrt((n+1)*(n+2))*(4*n+6)*kronecker(m,n+2) + sqrt((n + 1)*(n + 2)*(n + 3)*(n + 4))*kronecker(m,n+4))*mu + (n + 1/2)*kronecker(m,n);
    end
end

%eigenvalues and eigenvectors
[res2,res1] = eig(ham);
evals=diag(res1);
evecs=res2;
%energy eigenvalues
evals
%overlaps <psi_j|0> 
evecs(1,:)

%test that we recover the t = 0 result
x=(-40:40)/10;
y1=zeros(numel(x),1);
y2=zeros(numel(x),1);
for q = 1:numel(x)
        y1(q)=my_p(x(q),0,nn,evecs,evals);
        y2(q)=abs(my_psi(0,x(q)))^2;
end

plot(x,y1,'o')
hold on
plot(x,y2,'-')

%% 
%other values of t
x=(-30:30)/10;
res05=zeros(numel(x),1);
res1=zeros(numel(x),1);
res15=zeros(numel(x),1);
res2=zeros(numel(x),1);
for q = 1:numel(x)
        res05(q)=my_p(x(q),0.5,nn,evecs,evals);
        res1(q)=my_p(x(q),1.0,nn,evecs,evals);
        res15(q)=my_p(x(q),1.5,nn,evecs,evals);
        res2(q)=my_p(x(q),2.0,nn,evecs,evals);
end

plot(x,res05,'o')
hold on
plot(x,res1,'o')
hold on
plot(x,res15,'o')
hold on
plot(x,res2,'o')



%FUNCTION DEFINITIONS
function y=kronecker(a,b) 
    if a==b
        y=1;
    else
        y=0;
    end
end

%Harmonic oscillator wave functions in dimensionless variables
function y=my_psi(n,x)
    y = 1/sqrt(factorial(n)*2^n)*hermiteH(n,x/sqrt(2))*exp(-x^2/4)/(2*pi)^(1/4);
end

%dimensionless probability density
function y=my_p(z,tau,nn,vec,val)
    y=0.0;
    for n = 1:2:(nn+1)
        for k = 1:2:(nn+1)
            y=y+my_psi(k-1,z)*exp(-1i*tau*val(n))*vec(1,n)*vec(k,n);
        end
     end
    y=abs(y)^2;
end
%here we sum only over half the modes because the other overlaps are zero due to parity symmetry 


