QM4(b)

clc

clf

m=940e6

h=1973

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

n=2500

T=zeros(n-1,n-1);

for i = 1:n - 1

T(i,i)= - 2;

if i > 1

T(i, i - 1)= 1

end

if i < n - 1

T(i, i + 1)= 1 ;

end

end

xmin =0;

xmax =(30);

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

x=xmin+dx:dx:xmax

xp=(x-x0)./x

v=D*(exp(-2*a*xp)-exp(-a*xp))

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(1,2,1)

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

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

xlabel("Position(x1)")

ylabel("Unl_state")

a = gca();

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

subplot(1,2,2)

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;10, 0.27];