%% For 2 Eigenstates

close all; clear all; clc
a = 1
N = 100
x = linspace(0,a,N);
t = linspace(0,20*pi,200)
E2 = 12.9
E1 = 13.6

c = 1/sqrt(2);
Psi1 = sqrt(2/a)*sin(1*pi*x/a);
Psi2 = sqrt(2/a)*sin(2*pi*x/a);
Psinet = c*Psi1+c*Psi2;
subplot(3,1,1)
plot(x, Psi1,'-', x, Psi2,'--', x, Psinet,'*')
legend('Psi1(x,0)', 'Psi2(x,0)', 'Psi Net(x,0)','location','best')
xlabel('x')
ylabel('Psi')

subplot(3,1,2)
plot(x, Psi1.^2,'-', x, Psi2.^2,'--', x, (Psi1.^2+Psi2.^2),'*')
xlabel('x')
ylabel('|Psi(x,0)|^2')
legend('|Psi1(x,0)|^2', '|Psi2(x,0)|^2', '|Psi Net(x,0)|^2','location','best')

% Wave Function Oscillation
PsiT0 = c^2*(Psi1.^2 + Psi2.^2+2.*Psi1.*Psi2.*cos((E2-E1)*t(1)));
A = axis
subplot(3,1,3)
linehandle = plot(x,PsiT0, '--')
axis(A)
xlabel('x')
ylabel('|Psi(x,t)|^2')
Title('Constructive and Destructive Oscillations')
for i = 1:length(t)
    PsiT = c^2*(Psi1.^2 + Psi2.^2+2.*Psi1.*Psi2.*cos((E2-E1)*t(i)));
    set(linehandle, 'ydata',PsiT);
    drawnow
    pause(.05)
end

%% Now for N Eigenstates
close all; clear all; clc
%inputs
L = 1;
A = 500;
x = linspace(0,L,A);
t = linspace(0,20*pi,2000);
N = 2;
E0 = 13.6;
E = linspace(E0, 1, N);
c = 1/sqrt(N);

for i= 1:N
    PsiNet = zeros(size(x));
    Psi{i} = c*sqrt(2/L)*sin(i*pi*x/L);
end
clear i
for i = 2:N
    PsiNet = Psi{1} + Psi{i};
end

subplot(2,1,1);
plot(x, PsiNet,'--r','linewidth',2);
xlabel('x');
ylabel('Psi(x,0)');
title(sprintf('N=%d',N));


%now we construct time dependent ones
PsiT = cell(N, length(t));

clear i
for i = 1:length(t)
    for k = 1:N
        PsiT{k,i}=Psi{k}.*exp(-j.*E(k).*t(i));%Psi is in Rows, time is in columns
    end
end
clear i
clear k


PsiTNet = cell(1,length(PsiT));

for k = 1:length(PsiT);
    for i = 2:N
        PsiTNet{k} = PsiT{1,k}+PsiT{i,k}; %This is a row vector which contains
        %the wave functions, with time information allong the column.
    end
end
%now I have to make the Probability distribution
PsiTNetSquared = cell(size(PsiTNet));
for i = 1:length(PsiTNet);
    PsiTNetSquared{i} = PsiTNet{i}.*conj(PsiTNet{i});
end

%now I plot this business in an animation to see if it worked
hold on
subplot(2,1,1)
plot(x,PsiTNetSquared{1}, 'k', 'linewidth', 2)
ylabel('')
legend('Psi(x,0)','|Psi(x,0)|^2')


subplot(2,1,2);
linehandle = plot(x, PsiTNetSquared{1},'k', 'linewidth',1.5);
xlabel('x');
ylabel('|Psi(x,t)|^2');
clear i;
A = axis;
axis(A) % i'm doing this so that the window doesnt move around
for i=1:length(PsiTNetSquared)
    set(linehandle, 'ydata', PsiTNetSquared{i});
    drawnow
    pause(.1)
    title(sprintf('N = %d, Time = %d', N, t(i)))
end