ps8 p2.  Compute the matrix element for the p-orbital (l=1) subset of V = λ ((x^2)-(y^2)), diagonalize it, and find the diagonalizing basis states.

ClearAll[b, braket, theta, phi, display, sopts, v, vm, diagYdisplay, diagY, energies, states]

b = {SphericalHarmonicY[1,1,theta,phi],
SphericalHarmonicY[1,0,theta,phi],
SphericalHarmonicY[1,-1,theta,phi]} ;
y = {"(Subscript[Y, 1, 1])(θ, ϕ)", "(Subscript[Y, 1, 0])(θ, ϕ)","(Subscript[Y, 1, -1])(θ, ϕ)" } ;
v = Sin[theta]^2 (Cos[phi]^2 - Sin[phi]^2) ;
display = {theta -> θ, phi -> ϕ} ;

braket[f_, g_] := Integrate[ Conjugate[f] g Sin[theta],{theta, 0, Pi}, {phi, 0, 2 Pi}] ;

row[n_] :=braket[b[[n]],v b[[#]]] &/@ Range[3] ;
vm = row[#] &/@Range[3] ;

{energies, states} = Eigensystem[ vm ] ;
states = (#/Norm[#]) &/@ states ;

diagYdisplay =(#.y ) &/@states //Simplify ;
diagY =(#.b ) &/@states // TrigToExp // FullSimplify ;

sopts = {14, "Text"} ;
Style[{#1," = ", #2 } &@@@{
{"m = 1,0,-1 ; (Subscript[Y, 1, m])", b /. display// Column // TraditionalForm  },
{"<n'1m'|V|n1m>/(∫_0^∞)(Subscript[R, n]^2) (r) (r^4) λ\[DifferentialD]r",vm // MatrixForm},
{"Δ\!\(\*TagBox[
RowBox[{"/", 
RowBox[{SubsuperscriptBox["∫", "0", "∞"], 
RowBox[{SubsuperscriptBox["R", "n", "2"], " ", 
RowBox[{"(", "r", ")"}], " ", SuperscriptBox["r", "4"], " ", "λ", 
RowBox[{"\[DifferentialD]", "r"}]}]}]}],
HoldForm]\)",energies // Column },
{"diagonalizing basis", {(# // MatrixForm)  ,","} &/@states // Flatten // Row },
{"<θ, ϕ|ψ>",diagYdisplay // Column},
{"<θ, ϕ|ψ>",(# /. display// TraditionalForm) &/@ diagY // Column}
} // Grid ,sopts]