Spin three halves operator question (Sakurai 3.33).  Implement operators for (SubPlus[S]), (SubMinus[S]), and (Subscript[S, z]) that act on kets |±3/2>,|±1/2>.  These operators act on only a single basis element, but are used to construct the matrix representations for these operators (which are more general).  Use those matrices to compute the representation of the Hamiltonian for the problem and compute its energy eigenvalues.

ClearAll[a,b, hbar, display, basisVectors, mValues, m, vec, raiseOrLower, raise, lower, splus, sminus, sz, szOp, sSquared, h, sopts]
On[Assert] ;
$Assumptions=(And@@(#>0&)/@{hbar, a,b});

display := {hbar -> ℏ, a-> "A", b->"B"} ;
basisVectors = IdentityMatrix[4] ;

(*values for spin 3/2*)
mValues := {3/2,1/2,-1/2,-3/2} ;


(*Map the values m, in|m> to the standard basis for (R^4)*)
vec[3/2] := basisVectors[[1,All]] ;
vec[1/2] := basisVectors[[2,All]] ;
vec[-1/2] := basisVectors[[3,All]] ;
vec[-3/2] := basisVectors[[4,All]] ;

(* Any invalid value is mapped to a zero vector (i.e. (SubPlus[S])|3/2>, or (SubMinus[S])|-3/2>*)
vec[n_] := {0,0,0,0} ;

(* For a basis vector |m> (not a linear combination) return the numeric value m *)
m[v_] := v.DiagonalMatrix[mValues] // Tr ;

(*s = 1 for raise*)
raiseOrLower[v_, s_Integer] := Module[{mv},
Assert[s == 1 || s == -1] ;
mv = m[v] ;
hbar Sqrt[(3/2 -s mv)(5/2 + s mv)] vec[mv+ s ]
] ;

raise[v_] := raiseOrLower[v,1] ;
lower[v_] := raiseOrLower[v,-1] ;
szOp[v_] := hbar m[v] v ;

(* Form the matrix for an operator that acts on vec[m] values. *)
operatorToMatrix[o_] := (o[vec[#]] &/@mValues)// Transpose ;

(* Matrixes for (SubPlus[S]), (SubMinus[S]), and (Subscript[S, z]), (S^2) *)
splus = operatorToMatrix[raise] ;
sminus = operatorToMatrix[lower] ;
sz := operatorToMatrix[szOp] ;
sSquared := (3/2)(3/2+1) hbar^2 IdentityMatrix[4] ;

(* The Hamiltonian for the problem. *)
h = a ( 3 sz.sz - sSquared) + b ( splus.splus  + sminus.sminus) ;

(* Display the Hamiltonian and the energy eigenvalues. *)
sopts = {14, "Text"} ;

Style[{#1," = ", #2  /. display } &@@@{
{"(SubPlus[S])",splus // MatrixForm },
{"(SubMinus[S])",sminus  // MatrixForm },
{"(Subscript[S, x])",(splus +  sminus)/(2)  // MatrixForm },
{"(Subscript[S, y])",(splus - sminus)/(2I) // MatrixForm },
{"(Subscript[S, z])",sz // MatrixForm },
{"H",h  // MatrixForm },
{"eigenvalues",Eigenvalues[h] // Column}
} // Grid ,sopts]

Scratch expressions.

m[vec[3/2]]
(vec[#] // MatrixForm)&/@ mValues

(raise[vec[#]] // MatrixForm)&/@ mValues
(lower[vec[#]] // MatrixForm)&/@ mValues
splus // MatrixForm
sminus // MatrixForm

(splus . vec[#] // MatrixForm)&/@ mValues
(sminus . vec[#] // MatrixForm)&/@ mValues

splus.splus // MatrixForm
sminus.sminus // MatrixForm
sz.sz // MatrixForm
sSquared // MatrixForm