Clear[ norm, expectation, phi, r, z, f, θ]

phi[r_] := (z/(Subscript[a, 0 ]))^(3/2) (2/Sqrt[4 Pi]) E^(-z r/(Subscript[a, 0]) )       (* l = 0, S state ; n=1: 1S state  *)
expectation[f_] := 2 Pi Integrate[ Integrate[f Sin[θ], {θ, 0, Pi}] r^2 (phi[r])^2, {r, 0, Infinity}]

(* This is for verification of 16.74 and 16.75 *)

expectationZsquared := expectation[ r^2 (Cos[θ])^2]
expectationRtimesZsquared := expectation[ r (Cos[θ])^2 r^2]
expectationRcubed := expectation[ r^3]
expectationRsquared := expectation[ r^2]
(*
expectationZsquared/expectationRsquared
expectationRtimesZsquared/expectationRcubed 
*)
(e^2 absE^2)(-m (Subscript[a, 0]) )/(ℏ^2)expectation[(r/2 + (Subscript[a, 0])) r^2 (Cos[θ])^2] /. {z -> 1, (m/ℏ^2) -> (1/(e^2 (Subscript[a, 0])))}

(* An attempt to verify 16.72 *)

Clear[f, r, θ]
f[r_, θ_] = -(m aNought/ℏ^2)(r/2 +aNought) r Cos[θ]
Clear[laplacian]
laplacian[g_] := (1/r^2) D[r^2 D[ g, r], r]  + (1/(r^2 Sin[θ])) D[ Sin[θ] D[ g, θ]] 
(Subscript[h, 0]) := (-ℏ^2/(2m))
SetAttributes[{m, aNought, ℏ},Constant]
FullSimplify[f[r, θ] (Subscript[h, 0 ])laplacian[ψ[r, θ]]  -  (Subscript[h, 0 ])laplacian[ψ[r, θ] f[r, θ]]   ]
laplacian[ (r/2 + aNought) r Cos[θ]] // FullSimplify

(1/(4 r))aNought ((-2 (2 aNought+3 r) Cos[θ]+(2 aNought+r) Sin[θ]) ψ[r,θ]-4 r (aNought+r) Cos[θ] Derivative[1,0][ψ][r,θ]) /. aNought -> (Subscript[a, 0])

(((Subscript[a, 0]) ((Sin[θ] (r+2 (Subscript[a, 0]))-2 Cos[θ] (3 r+2 (Subscript[a, 0]))) ψ[r,θ]-4 r Cos[θ] (r+(Subscript[a, 0])) Derivative[1,0][ψ][r,θ]))/(4 r)) //FullForm

expectation[ (r/2 + (Subscript[a, 0 ])) r Cos[θ]]