MATLAB Code: Stationary Schrodinger EQ, E Spec, Eigenvalues

[Baynie]

Hello everyone,

For weeks I have been struggling with this quantum mechanics homework involving writing a code to determine the energy spectrum and eigenvalues for the stationary Schrödinger equation for the harmonic oscillator. I can't find any resources anywhere. If anyone could help me get started, get my matrices and equations set up, or has worked a similar problem/written a similar code before, any help would be greatly appreciated! Thanks in advance!

Homework Statement

Homework Equations

Included in image above

The Attempt at a Solution

Her is the code I have written so far. I'm not sure if this is even close or on the right track
% Stationary Schrödinger Equation - QHO

clear
clc

hbar = 6.58E-16;
f = 400E-9;
w = 2*pi*f;
m = 1;
N = 101;
a = 0.1;
n = 1:N;
r = n*a;
l = 1;
x = a;
mwhbar = m*w*hbar;
y = r;

% Operators (Matrices)
T = diag(-1*ones(1,N-2),2) + diag(2*ones(1,N-1),1) + diag(-1*ones(1,N),0);
K = (1/(2*a^2)) * T; % Kinetic Energy Matrix

Veff = -(1./r) + l*(l + 1)./(2*(r.^2));
V = diag(Veff);
U = (1/(2*a^2)) * V; % Potential Energy Matrix

% Equations
H = -(1/2*a^2)*(eigen_f(n+1) - 2 * eigen_f(n) * eigen_f(n-1));

 

[Baynie]

UPDATE:

Here is an update on the code I've been working on. It is probably closer to where I'm supposed to be headed. I just have problems when it comes to calculating V and Veff clear
clc

N = 101;
a = .01;
n = 1:101;
r(1:N) = n.*a;
l = input('Enter angular quantum number l: ');
T = zeros(N);
V = zeros(N);
Veff = zeros(N);for i = 1:N
for j = 1:N
if i == j
T(i,j) = 2/(2*a^2);
elseif i == j-1 || j == i-1
T(i,j) = -1/(2*a^2);
end
end
endfor i = 1:N
Veff(i,i) = (1/r) + (l*(l+1))./(2.*r.^2);
V(i,i) = V(i,i).*Veff;
end

H = T + V;

[Psi,Energy] = eig(H);

E = diag(Energy);

fprintf('Lowest Eigenvalues:\n');
disp(E(1:3))Psi_0 = Psi(:,1);
Psi_0_Squared = Psi_0.^2;

Psi_1 = Psi(:,2);
Psi_1_Squared = Psi_1.^2;Range_a = a:a:N*a;Integral_0 = trapz(Range_a,Psi_0_Squared);
Integral_1 = trapz(Range_a,Psi_1_Squared);Normalized_Psi_0 = 1/sqrt(Integral_0)*Psi_0;
Normalized_Psi_1 = 1/sqrt(Integral_1)*Psi_1;

subplot(3,1,1);
plot(Range_a,Normalized_Psi_0,'r');
title('Wave Function: Ground State Hydrogen Atom');
xlabel('x');
ylabel('Psi_0(x)');subplot(3,1,2);
plot(Range_a,Normalized_Psi_1,'r');
title('Wave Function: First Excited State Hydrogen Atom');
xlabel('x');
ylabel('Psi_1(x)');subplot(3,1,3);
plot(Range_a,Normalized_Psi_1,'r');
title('Wave Function: First Excited State Hydrogen Atom');
xlabel('x');
ylabel('Psi_2(x)');