Variational (Subscript[H, 2]^(+))

taking ψ=((Subscript[P, 1,0])((Subscript[r, A]),Z)+(Subscript[P, 1,0])((Subscript[r, B])),Z)(e^(-(z^2)/(a^2)))

(Subscript[P, {n_,l_}])[r_,Z_]:=(Sqrt[(((n-l-1)!Z )/((n^2)((n+l)!)))])(((((2 Z r)/n)))^(l+1))(E^(-((Z r)/n)))LaguerreL[n-l-1,2l+1,((2Z r)/n)]
(Subscript[Rnl, {n_,l_}])[r_,Z_]:=(((Subscript[P, {n,l}])[r,Z])/r)

((⟨ψ|H|ψ⟩)/(⟨ψ|ψ⟩))

energy[R_,Z_,a_]:=With[{coord={rA->(Sqrt[(((z+R/2))^2)+(ρ^2)]),rB-> (Sqrt[(((z-R/2))^2)+(ρ^2)])}},
psi=((Subscript[Rnl, {1,0}])[rA,Z]+(Subscript[Rnl, {1,0}])[rB,Z])(E^(-(z^2)/(a^2)))//.coord;
g_~hh~f_:=g (((-1)/2)((1/ρ)(∂_ρ)(ρ (∂_ρ)f)+(∂_(z,z))(f ))+((1/R)-(1/rA)-(1/rB))f )//.coord;
NIntegrate[ρ(psi~hh~psi),{ρ,0,∞},{z,-∞,∞}]/NIntegrate[ρ psi psi,{ρ,0,∞},{z,-∞,∞}]]

find minimum; first vary a:

Table[{a,energy[2,1,a]},{a,8}]

now vary Z at a=4

Table[{Z,energy[2,Z,4]},{Z,.8,1.5,.1}]

redo a

Table[{a,energy[2,1.2,a]},{a,2,5,.2}]

energy[2,1.2,3.8]

Plot[psi/.{ρ->0},{z,-5,5},AspectRatio->.3]