%NANO 705 - Spring 2014
%Homework 2 - Problem 3

clear all

%Constant
hc=1240;%eV.nm
mc2=511e3;%eV

%Parameters
L=10;%nm
Nx=10000;%divisions in x
periodic_BCs=0;%1: true, 0: false

%Ng=5;Vgmin=0.1;Vgmax=0.5;
xc=linspace(0.5*L/Nx,L-0.5*L/Nx,Nx);
Psi=[1:Nx];%=linspace(0,1.0,Nd);
fac=(hc/2/pi)^2/2/mc2;

for ix=1:Nx
    Psi(ix)=(2/L)^0.5*sin(pi*xc(ix)/L);%eigen function
    %Psi(ix)=(30/L)^0.5*(xc(ix)/L)*(1-xc(ix)/L);%trial function
end
dx=L/Nx;

disp('------------------------------');
format long

%Check normalization
Itot=0;%integral
for ix=1:Nx
    Itot=Itot+Psi(ix)^2*dx;
end
Itot

%Find energy expectation value
Etot=0;%energy (eV)
for ix=1:Nx
    im=ix-1;
    ip=ix+1;
    if(periodic_BCs)
        if(im<1)
            im=Nx;
        end
        Psim=Psi(im);
        
        if(ip>Nx)
            ip=1;%
        end
        Psip=Psi(ip);
        
    else
        if(im<1)
          Psim=0;
        else
          Psim=Psi(im);
        end
 
        if(ip>Nx)
          Psip=0;
        else
           Psip=Psi(ip);
        end
        
    end
    d2Psi_dx2=(Psip-2*Psi(ix)+Psim)/(dx^2);
    HPsi=-fac*d2Psi_dx2;
    Etot=Etot+Psi(ix)*HPsi*dx;
end
Etot

hold on;
p=plot(xc,Psi);
set(p,'linewidth',[4.0])
set(gca,'Fontsize',[12])
xlabel('x (nm)')
ylabel('Psi')
%legend('eps-mu','eps(0)-mu')
grid on