Clear[a,  h, r, norm, e, eu, ψ]
h[r_] := (1/(Pi a^2))^(3/4)E^(-r^2/(2a^2))
norm = 4 Pi Integrate[ r^2  h[r] h[r] , {r, 0, Infinity}]
e := 4 Pi  Integrate[ r^2 ( h[r] (-ℏ^2/(2m)) (h''[r] + 2 h'[r]/r) + h[r] h[r] (1/2) m ω^2 r^2), {r, 0, Infinity}]
eu := e /. a -> Sqrt[u] 
eu = FullSimplify[eu, Re[u] > 0]
Solve[D[eu, u] == 0, u]
(* subst the positive root back in the variational energy *)
FullSimplify[eu /. u -> ℏ/(m ω)]
ψ[r_] =h[r] /. a -> Sqrt[(ℏ/(m ω))] (* extremum wavefunction *)
FullSimplify[(-ℏ^2/(2m)) (ψ''[r] + 2 ψ'[r]/r) +  ψ[r] (1/2) m ω^2 r^2]/ψ[r] (* verify that we get the same as the subst back into the variational energy *)