a=((ℏ^2)/(me k (qe^2)))
En[n_]:=-((k (qe^2))/(2 a))(1/(n^2))
β[n_]:=(Sqrt[-8 me ((En[n])/(ℏ^2))])
Psi[r_,theta_,phi_,n_,l_,m_]:=(Sqrt[((((2/(n a))))^3)(((n-l-1)!)/(2n (n+l)!))])Exp[-β[n] r/2](((β[n] r))^l)LaguerreL[n-l-1, 2l+1,β[n] r]SphericalHarmonicY[l,m,theta,phi];

Psi[r,θ,ϕ,1,0,0]
Psi[r,θ,ϕ,2,1,0]

(((E^(-(1/2) r (Sqrt[(((k^2) (me^2) (qe^4))/(ℏ^4))]))) r (Sqrt[(((k^3) (me^3) ((qθe)^6))/(ℏ^6))]) (Sqrt[(((k^2) (me^2) (qe^4))/(ℏ^4))]) Cos[θ])/(4 (Sqrt[2 π])))

Integrate[Cos[θ]Sin[θ],{θ,0,Pi}]

Integrate[(((E^(-2r (Sqrt[(((k^2) (me^2) (qe^4))/(ℏ^4))]))) (((k^3) (me^3) (qe^6))/(ℏ^6)))/π)(r^2),{r,0,∞},Assumptions->{k>0,me>0,qe>0,ℏ>0}]

Integrate[Conjugate[Psi[r,θ,ϕ,2,1,-1]]Psi[r,θ,ϕ,2,1,1](r^2)Sin[θ],{θ,0,π},{ϕ,0,2 π},{r,0,∞},Assumptions->{k>0,me>0,qe>0,ℏ>0}]

Simplify[(-((ℏ^2)/(2 me ))(D[Psi[r,θ,ϕ,1,0,0],{r,2}]+(2/r)D[Psi[r,θ,ϕ,1,0,0],{r,1}])-((k (qe^2))/r)Psi[r,θ,ϕ,1,0,0])/Psi[r,θ,ϕ,1,0,0],Assumptions->{k>0,me>0,qe>0,ℏ>0}]

En[1]