QM4(a)

clc

clf

m=940

e=3.795

h=197.3

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

n=1610

k=100

b_all=[0 10 30];

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

b=b_all(i)

v=1/2*(k*x^2)+ 1/3*(b*x^3);

PE= diag (v)

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

[evect,evalue] = spec(H)

Energy= diag(evalue)

disp("a="+string(a_f))

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')

//ax = gca();

//ax.data_bounds = [0, -0.15;20, 0.05];

title("a="+ string(a_f))

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

ylabel("Unl_state")

subplot(2,3,i+3)

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

//ax = gca();

//ax.data_bounds = [0, 0;20, 0.25];

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

xlabel("Position(A)")

ylabel("Probability Density")

end