Sakurai pr. 5.17 (a): H = A (L^2)+ B (Subscript[L, z])+ C (Subscript[L, y])

ClearAll[hamiltonian, h1, h2,h3, a, b, c, hbar, ℏ, sortbybFilter, sortByB ,displayEigenvalues] ;

$Assumptions=(And@@(#>0&)/@{hbar, a, b,c});

display = {hbar -> ℏ, a -> "A", b->"B", c-> "C"} ;

hamiltonian[l_Integer, sparse_ : False ] := Module[{r,h, d,lx},
r = 2 l ;
d = a l (l+1) hbar^2 + b hbar( -l + # -1) &/@ Range[r + 1] ;
lx = hbar Sqrt[#(r-# +1)] &/@ Range[r] ;

If [ sparse,
h = SparseArray[{
Band[{1,1}]->d,
Band[{2,1}]-> c  lx/I,
Band[{1,2}]-> -c  lx/I},{r+1,r+1}]
,
h = DiagonalMatrix[ d ] 
+ DiagonalMatrix[ c  lx/I, -1 ] 
- DiagonalMatrix[ c  lx/I, 1 ] ;
] ;

h
] ;

sortbybFilter[v_] := Module[{r},
r = (v /. {a->0, c->0, b->1})/hbar ;
(*r // Print ;*)
r] ;
sortByB[s_] := Sort[s, sortbybFilter[#1] > sortbybFilter[#2] &]  ;

displayEigenvalues[l_] := Module[{hl,e,s},
hl = hamiltonian[l] ;
e= Eigenvalues[hl] ;

s = ((Series[#, {c,0,2}] // Simplify // Normal // Expand)&/@ e );

Grid[{e // sortByB,s // sortByB} // Transpose] /. display
] ;

(hamiltonian[5] // MatrixForm) /. display
displayEigenvalues[5]

delta[l_Integer] :=-(c^2 hbar/(2b))(l^2-(-l+#) + l) &/@ (Range[2l+1]-1) ;

(delta[5] b/c^2/hbar // Simplify)// N