2D harmonic oscillator in sylondric cordinates code

Question

I'm trying to solve the Hamiltonian matrix for a 2D isotropic harmonic oscillator in cylindrical coordinates using MATLAB, but I'm running into some issues. Below is my MATLAB code and a description of the problem I'm facing. Can someone help me how I add the term 0.5*((l^2)/(rho^2)) at the diagonal with l angular momentum.

code part 2

% Berekening van gewichtsfactoren mu
mu= sqrt((1 / 2))*sqrt((rho1).^2-(rho2).^2);
mu_i_plus_1= mu(2:end);
mu_i_min_1= [mu(1)];
mu_i_min_1= [mu_i_min_1, mu(1:end)];
% Stapgrootte h = rho(2) - rho(1);
% Berekening van de diagonale en niet-diagonale elementen van de Hamiltoniaan
H_diag = (-rho1-rho2)./(h.*(mu).^2) - (l.*2 ./ (rho(1:end-1).^2)) - (rho(1:end-1).^2);
H_off = (rho1(1:end-1))./(mu(1:end-1).*mu_i_plus_1.*h);
H_0 = -0.5*(diag(H_diag)+ diag(H_off(1: length(H_diag)-1),1) + diag(H_off(1: length(H_diag)-1),-1));