%Numerisk losning av TUSL for partikkel i boks
%
clear all;
%m=masse (SI)
hbar=1.05E-34;
m=9.11E-31;  %elektron
V0=0.00;
%N = antall diskrete posisjonsverdier i boksen
N=100;
%dz = skrittlengde; boksbredde = N*dz
dz=0.040E-9;
z=zeros(1,N);
V=zeros(1,N);
for i=1:N
    z(i)=dz*i*1E9;   %tabell med posisjonsverdier i nm
    V(i)=V0;         %V=0 i boksen, uendelig utenfor
end
% Diagonale elementer i Hamiltonmatrisen
d = (1/dz^2)*hbar^2/m + V;
% Ikke-diagonale elementer:
e = (-1/(2*dz^2))*hbar^2/m;
H = zeros(N,N);
% Lager den tridiagonale H, først diagonalen
for i=1:N
    for j=1:N
        if i==j
            H(i,j)=d(i);
        end
        if abs(i-j)==1
            H(i,j)=e;
        end
    end
end
%Se Matlab Help for dokumentasjon av eig(H)
[psi,w] = eig(H);
%Lagrer egenverdiene i tabellen evalues, enhet eV:
evalues = diag(w)/1.60E-19;
%Skriver ut i kommandovinduet de 4 laveste energiegenverdiene i eV:
evalues(1) 
evalues(2) 
evalues(3) 
evalues(4) 
% Plotter psi1,..,psi4
psi1=psi(:,1);
psi2=psi(:,2);
psi3=psi(:,3);
psi4=psi(:,4);
%Plotter disse 4 som funksjon av z.
figure(1);
plot(z,psi1,z,psi2,z,psi3,z,psi4);
title('Particle in 1D box');
xlabel('z (nm)');
ylabel('\psi');
figure(2);
hold on;
plot([0 z(N)],[evalues(1) evalues(1)]);
plot([0 z(N)],[evalues(2) evalues(2)]);
plot([0 z(N)],[evalues(3) evalues(3)]);
plot([0 z(N)],[evalues(4) evalues(4)]);
plot([0 z(N)],[evalues(5) evalues(5)]);
title('Energy levels');
ylabel('E (eV)');
hold off;
figure(3);
plot(z,abs(psi1.^2),z,abs(psi2.^2),z,abs(psi3.^2),z,abs(psi4.^2));
title('Probability densities');
xlabel('z (nm)');
ylabel('|\psi|^2');