Hamiltonian for the (Subscript[H, 2]^(+)) molecule

H=((-1)/2)(∇^2)-(1/(Subscript[r, A]))-(1/(Subscript[r, B]))+(1/R)

Try (unnormalized) wavefunctions |n}=((|(Subscript[A, n])⟩+|(Subscript[B, n])⟩)/(Sqrt[2])) where |(Subscript[A, n])⟩=(((Subscript[P, nl])((Subscript[r, A])))/(Subscript[r, A]))(Subscript[Y, lm])((Subscript[Ω, A])) is a normalized hydrogenic wavefunction centered on nucleus A, (Subscript[B, n]) likewise centered on B. 

The Hamiltonian can be written in a convenient, if unsymmetrical form, as H=(Subscript[H, A]^0)+(1/R)-(1/(Subscript[r, B])), where (Subscript[H, A]^0)
is the hydrogen atom Hamiltonian, centered at atom A.  Thus (Subscript[H, A]^0)|(Subscript[A, n,l])⟩=(Subscript[E, n]^0)|(Subscript[A, n,l])⟩.

 The matrix elements of H are:
 {n|H|m}=(1/2)(⟨(Subscript[A, n])|H|(Subscript[A, m])⟩+⟨(Subscript[B, n])|H|(Subscript[B, m])⟩+⟨(Subscript[B, n])|H|(Subscript[A, m])⟩+⟨(Subscript[A, n])|H|(Subscript[B, m])⟩).  Symmetry requires ⟨(Subscript[A, n])|H|(Subscript[A, m])⟩=⟨(Subscript[B, n])|H|(Subscript[B, m])⟩ and ⟨(Subscript[B, n])|H|(Subscript[A, m])⟩=⟨(Subscript[A, n])|H|(Subscript[B, m])⟩, so
 {n|H|m}=⟨(Subscript[A, n])|H|(Subscript[A, m])⟩+⟨(Subscript[A, n])|H|(Subscript[B, m])⟩.
 
 Evaluating these matrix elements then gives
{n|H|m}=((Subscript[E, n]^0)+1/R)((Subscript[δ, n,m])+⟨(Subscript[A, n])|(Subscript[B, m])⟩)-⟨(Subscript[A, n])|1/(Subscript[r, B])|(Subscript[A, m])⟩-⟨(Subscript[A, n])|1/(Subscript[r, B])|(Subscript[B, m])⟩

Now, to deal with the non-orthgonality and un-normalized nature of the basis states.  We expand |ψ⟩=(∑_m)(Subscript[c, m])|m} and substitute into the Schrodinger equation:
{n|H|m}(Subscript[c, m])=E{n|m}(Subscript[c, m]) using the summation convention.  In matrix terms, this is
(Subscript[H, nm])(Subscript[c, nm])=E (Subscript[O, nm])(Subscript[c, m]), where (Subscript[O, nm])={n|m}=(Subscript[δ, n,m])+⟨(Subscript[A, n])|(Subscript[B, m])⟩ is a matrix of overlap integrals.  Then the eigenvalue equation  becomes
(O^(-1))H|(Subscript[u, p]))=E |(Subscript[u, p])).  Since (O^(-1))H is not Hermitian, its eigenstates |(Subscript[u, p])) are not orthogonal in the usual sense.  Instead, we define their dual vectors ((Subscript[v, q])| such that ((Subscript[v, q])|(Subscript[u, p]))=(Subscript[δ, p,q]).  We can easily find the dual vectors by taking the inverse of the (transpose) of the list of eigenvectors produced by Mathematica.

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

st={{1,0},{2,0},{2,1},{3,0},{3,1},{3,2},{4,0},{4,1},{4,2},{4,3}};
Rs={1.5,2.,2.5,3,4,6,9}

Table[st={{1,0}};
int[f_]:=NIntegrate[(ρ/2) (f//.{rA->(Sqrt[(((z+R/2))^2)+(ρ^2)]),rB-> (Sqrt[(((z-R/2))^2)+(ρ^2)])}),{ρ,0,∞},{z,-∞,∞}];
Do[𝕐A[l]=(Sqrt[4π])SphericalHarmonicY[l,0,th,0]/.Cos[th]->(((z+R/2))/rA);𝕐B[l]=(((-1))^l)(Sqrt[4π])SphericalHarmonicY[l,0,th,0]/.Cos[th]->(((z-R/2))/rB);
,{l,0,5}];
(Subscript[B, {n_,l_}]):=(Subscript[Rnl, {n,l}])[rB]𝕐B[l];
(Subscript[A, {n_,l_}]):=(Subscript[Rnl, {n,l}])[rA]𝕐A[l];
Olap=Table[δ[m,n]+int[(Subscript[A, st[[m]]])(Subscript[B, st[[n]]])],{m,Length[st]},{n,Length[st]}];
h=Table[(((-1)/(2 ((st[[n,1]])^2)))+(1/R)) Olap[[m,n]]-int[(((Subscript[A, st[[m]]])(Subscript[A, st[[n]]]))/rB)]-int[(((Subscript[B, st[[m]]])(Subscript[A, st[[n]]]))/rB)]
,{m,Length[st]},{n,Length[st]}];
{eval,evec}=Eigensystem[Inverse[Olap].h];in=Ordering[eval];eval=eval[[in]];evec=evec[[in]];Print[{R,eval}];
{R,Min[eval]}
,{R,Rs}]

ListPlot[%,Joined->True]

(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)

Table[With[{R=2.2},st={{1,0}};
int[f_]:=NIntegrate[(ρ/2) (f//.{rA->(Sqrt[(((z+R/2))^2)+(ρ^2)]),rB-> (Sqrt[(((z-R/2))^2)+(ρ^2)])}),{ρ,0,∞},{z,-∞,∞}];
Do[𝕐A[l]=(Sqrt[4π])SphericalHarmonicY[l,0,th,0]/.Cos[th]->(((z+R/2))/rA);𝕐B[l]=(((-1))^l)(Sqrt[4π])SphericalHarmonicY[l,0,th,0]/.Cos[th]->(((z-R/2))/rB);
,{l,0,5}];
(Subscript[B, {n_,l_}]):=(Subscript[Rnl, {n,l}])[rB,Z]𝕐B[l];
(Subscript[A, {n_,l_}]):=(Subscript[Rnl, {n,l}])[rA,Z]𝕐A[l];
Olap=Table[δ[m,n]+int[(Subscript[A, st[[m]]])(Subscript[B, st[[n]]])],{m,Length[st]},{n,Length[st]}];
h=Table[((((Z(Z-2)))/(2 ((st[[n,1]])^2)))+(1/R)) Olap[[m,n]]-int[(((Subscript[A, st[[m]]])(Subscript[A, st[[n]]]))/rB)]-int[(((Subscript[B, st[[m]]])(Subscript[A, st[[n]]]))/rB)]
,{m,Length[st]},{n,Length[st]}];
{eval,evec}=Eigensystem[Inverse[Olap].h];in=Ordering[eval];eval=eval[[in]];evec=evec[[in]];Print[{Z,eval}]];
{Z,Min[eval]}
,{Z,0.5,1.5,.2}]

ListPlot[%,Joined->True]