// Submitted by Dr. Triranjita Srivastava. Assistant Professor, Physics Dept., Kalindi College, University of Delhi
//Operating system: Windows 
//SCILAB Ver: 5.5.2

//Objective: Application to solve differential equations for a bound system - Eigenvalue Problem

// Example:Let us find out the energy eigenvalues and corresponding wavefunction of a particle of mass 'M' trapped in infinite potential Well (potential V=0) of width 'L'
//We implement Finite Difference Method (FDM) to obtain the eigenvalues
// By using FDM the second order differential operator  is replaced by a trigonal matrix and the problem reduces to a simple eigenvlaue problem 


clc
clear
h_cut=1.05457*10^-34                //(Plancks constant/2pi) J-s
L=input("Enter the width of the potential well L (in m) = ")
M=input("Enter mass of particle M (in kg) = ")
n=250                               // Number of divisions for FDM
N=(2*n)+1
x1=0                               // Initial value of x-coordinate 
s=(L-x1)/N                         // Step size for implementing FDM
EV=6.242*10^18                      // joule to eV conversion
//******************************************************************************************
// Hamiltonion Matrix H=T+V;  T=Kinetic energy operator (-d^2/dx^2)*h_cut^2/2M);    V= 0 (for infinite potential well)
//******************************************************************************************
T=zeros(N-1,N-1)
for i=1:(N-1)
    x1=x1+s
    T(i,i)=2
    if (i<(N-1))
        T(i,i+1)=-1
        T(i+1,i)=-1
end
end

H=(T*h_cut^2*EV/(2*M*s^2))                        //Hamiltonion Matrix

//*************************************************************************
// Finding eigenvalues and corresponding wavefunctions
//*************************************************************************
eigenvalues=spec(H)
disp("The eigenvalues (eV) of three lowest states obtained by FDM are " )
disp(eigenvalues(1:3))
[U,z]=spec(H)

//*************************************************************************
// Plotiing of three lowest order wavefunctions
//*************************************************************************
x=linspace(s,L,N-1)                     // creating x-coordinates for potential well
xlabel('x-coordinate (10^-10 m)','fontsize',5)
ylabel('Wavefunction (a.u.)','fontsize',5)
title('Graph of Wavefunction for three lowest order mode','fontsize',5)
a=get("current_axes")      //get the handle of the newly created axes
a.font_size=2
t=get("hdl")                //get the handle of the newly created object
t.font_size=5    
plot(x*10^10,U(:,1)'./max(U(:,1)),'r','linewidth',3)
plot(x*10^10,U(:,2)'./max(U(:,2)),'b','linewidth',3)
plot(x*10^10,U(:,3)'./max(U(:,3)),'g','linewidth',3)
hl=legend(['Ground State';'I Excited State';'II Excited State'],5)
h1.font_size=2

//*************************************************************************
// Comparison of obtained eigenvalues with analytical solution
//*************************************************************************
disp("The eigenvalues (eV) of three lowest states obtained by analytical results are " )
for j=1:3
    E(j)=j^2*%pi^2*h_cut^2*EV/(2*M*L^2)
    disp (E(j))
end
