QM3(b)

clc

m=0.511D6

e=3.795

h=1973

s=-(h*h)/(2*m)

n=1610

a_all=[3 5 7];

T = diag(-2*ones(n-1,1)) + diag(ones(n-2,1),1) + diag(ones(n-2,1),-1);

xmin = 0;

xmax = 30;

dx= (xmax-xmin)/(n-1);

x=xmin+dx:dx:xmax

for i=1:3

a_f=a_all(i)

v=(-e^2./x)*(%e^(-r/a));

PE= diag (v)

H=(T* (s/dx^2)+PE)

[evect,evalue] = spec(H)

Energy= diag(evalue)

disp("n="+string(n))

disp (Energy(1:5))

x1=xmin:dx:xmax

Unl_gs=[0;evect(:,1)]

Unl_fs=[0;evect(:,2)]

subplot(2,3,i)

plot(x1, Unl_gs,'b--',x1, Unl_fs,'r')

legend("Ground state", "First.sate")

ylabel("Unl_state")

a = gca();

a.data_bounds = [0, -0.15;14, 0.1];

subplot(2,3,i+3)

plot(x1,4*%pi*Unl_gs.^2,'b--',x1,4*%pi*Unl_fs.^2,'r')

legend("Gs density","Fs density")

xlabel("Position(x1)")

ylabel("probability Density")

a=gca();

//a.data_bounds = [0, -0.05;10, 0.25];

xlabel("Position(x1)")

end