SciLab: Solving the s-wave Schrodinger equation for the ground state and the first excited state of the hydrogen atom.

 The following script uses the Finite Difference Method (FDM) to solve time-independent Schrodinger equation (TISE) to find out the eigenstates of H atom.

// Aim: Solve the s-wave Schrodinger equation for the ground state and the first excited

//state of the hydrogen atom.

//Here, m is the reduced mass of the electron. Obtain the energy eigenvalues and plot

// the corresponding wavefunctions. Remember that the ground state energy of the

//hydrogen atom is -13.6 eV. Take e = 3.795 (eVÅ)^1/2, ħc = 1973 (eVÅ) and m = 0.511x10^6 eV/c2.

close;

clear;

clc;

// declaring constant values

hbarc=1973;//Plancks constant h divided by 2*(pi) called as hbar=(h/2*pi). This factor when multiplied by speed of light c then hbar*c in the units of (eVÅ) comes out to be 1973;

mcsq=0.511*10^(6);//This is mass of electron*c ^2 we call it mcsq in units of (eV);

e = 3.795; // (eVÅ)^1/2

// We hereby input the 'r' values so as to obtain the potential V(r) as a function of 'r')

r_min=0 // in units of angstrom

r_max=18 // in units of angstrom

N =  input("Input the number of intervals (should be around 500 to 750 for good computation)")

s = (r_max-r_min)/N; //step size

factor1=-(hbarc^2)/(2*mcsq*s^2);  // this factor is (hbar^2*c^2/2m*c^2) divided by s^2 //k=(hbar_c*hbar_c)/(2*m)

// Kinetic energy matrix (Using central difference formula)

T=zeros(N-1,N-1)

for i=1:N-1

    T(i,i)=-2

    if i<N-1 then

        T(i,i+1)=1

        T(i+1,i)=1

     end

end

T_matrix = factor1*T;    // Kinetic Energy Matrix final (scaling done)

//Potential energy matrix

V_matrix=zeros(N-1,N-1)

for i=1:N-1

    r(i)=r_min+i*s

    V_matrix(i,i)=-(e*e)/r(i);

end

// Hamiltonian matrix

H_matrix=T_matrix+V_matrix

// energy eigenvalue and eigenstates

[u,eigen]=spec(H_matrix);

// displaying of the ground and first excited state energies

disp("Ground state energy (1S orbital) for hydrogen atom is (in eV) : ")

disp(eigen(1,1))

disp("First excited state (2S orbital) energy  for hydrogen atom is (in eV): ")

disp(eigen(2,2))

disp("Second excited state (3S orbital) energy  for hydrogen atom is (in eV): ")

disp(eigen(3,3))

disp("Third excited state (4S orbital) energy  for hydrogen atom is (in eV): ")

disp(eigen(4,4))

rmatrix = [0;r;18]; //including the first and last point at which wavefunction is zero. 

//Displaying the first four energy eigen values

figure(1);

//scf()

for n = 1:1:4

plot(rmatrix, eigen(n,n)*ones(N+1,1), 'linewidth',n)

hl=legend(['n = 1';'n= 2';'n=3';'n=4']);

title('Energy Eigen values for Hydrogen atom','fontsize',3)

xlabel('r (Angstrom)','fontsize',3)

ylabel('Eigen Value','fontsize',3)

end

//Plotting the Probability |Psi^2| as a function of r

figure(2);

//scf()

R_wave_1s=u(:,1)./r;   //Radial wavefuction

R_wave_1s_final=[0;R_wave_1s;0];  //including the first and last point at which wavefunction is zero. 

// plot of probability function 

P_wave_1s= R_wave_1s_final.*R_wave_1s_final; 

plot(rmatrix,P_wave_1s, 'linewidth',3)

title('Plot of Probability function for 1s orbital', 'fontsize',3)

xlabel('r (Angstrom)','fontsize',3)

ylabel('Probability','fontsize',3)

figure(3);

//scf()

R_wave_2s=u(:,2)./r;

R_wave_2s_final=[0;R_wave_2s;0];   //including the first and last point at which wavefunction is zero. 

// plot of probability function 

P_wave_2s= R_wave_2s_final.*R_wave_2s_final; 

plot(rmatrix,P_wave_2s, 'linewidth',3)

title('Plot of Probability function for 2s orbital','fontsize',3)

xlabel('r (Angstrom)','fontsize',3)

ylabel('Probability','fontsize',3)

//Plotting the Probability density |Psi^2|*r^2 as a function of r

figure(4);

//scf() 

plot(rmatrix,[P_wave_1s, P_wave_2s], 'linewidth',3)

title('Plot of Probability function for 1s and 2s orbital','fontsize',3)

xlabel('r (Angstrom)','fontsize',3)

ylabel('Probability','fontsize',3)

//Plotting the Probability of electron |Psi^2 r^2| as a function of r

figure(5);

//scf()

fp = 4*3.14

Rr_wave_1s=u(:,1).*fp;   //Radial wavefuction

Rr_wave_1s_final=[0;Rr_wave_1s;0];  //including the first and last point at which wavefunction is zero. 

// plot of probability function 

Pr_wave_1s= Rr_wave_1s_final.*Rr_wave_1s_final; 

plot(rmatrix,Pr_wave_1s, 'linewidth',3)

title('Plot of Probability Probability of finding the e for 1s orbital', 'fontsize',3)

xlabel('r (Angstrom)','fontsize',3)

ylabel('Probability','fontsize',3)

figure(6);

//scf()

Rr_wave_2s=u(:,2).*fp;

Rr_wave_2s_final=[0;Rr_wave_2s;0];   //including the first and last point at which wavefunction is zero. 

// plot of probability function 

Pr_wave_2s= Rr_wave_2s_final.*Rr_wave_2s_final; 

plot(rmatrix,Pr_wave_2s, 'linewidth',3)

title('Plot of Probability of finding the e for 2s orbital','fontsize',3)

xlabel('r (Angstrom)','fontsize',3)

ylabel('Probability','fontsize',3)

figure(7);

//scf()

Rr_wave_3s=u(:,3).*fp;

Rr_wave_3s_final=[0;Rr_wave_3s;0];   //including the first and last point at which wavefunction is zero. 

// plot of probability function 

Pr_wave_3s= Rr_wave_3s_final.*Rr_wave_3s_final; 

plot(rmatrix,Pr_wave_3s, 'linewidth',3)

title('Plot of Probability of finding the e for 3s orbital','fontsize',3)

xlabel('r (Angstrom)','fontsize',3)

ylabel('Probability','fontsize',3)