e = 1.602176487 10^(-19);
m0 = 9.10938215 10^(-31);
e0 = 8.85418781762 10^(-12);
hbar = 1.054571628 10^(-34);

a0 = ((4 Pi e0 (hbar^2))/(m0 (e^2)));

e1 = -((m0 (e^4))/(32 (Pi^2) (e0^2) (hbar^2)));

s[r_] := (1 + (r/a0) + ((r^2)/(3 (a0^2)))) Exp[-r/a0];

haa[r_] := e1 + ((e^2)/(4 Pi e0 a0))(1 + (a0/r))Exp[-2 r/a0]

hab[r_] := e1 s[r] + ((e^2)/(4 Pi e0 a0 )) ((s[r] a0)/r) - ((e^2)/(4 Pi e0 a0 ))(1 + (r/a0)) Exp[-r/a0]

eLCAO[r_] := ((haa[r] + hab[r])/(1 + s[r]));

Plot[eLCAO[r a0]/e, {r, 0,12}, PlotPoints -> {100}, PlotRange -> {{0, 12},{0, -16}}, AxesLabel->{"R/(Subscript[a, 0])","Energy/eV"}]

min = FindMinimum[eLCAO[r a0]/e, r]

rMin = 2.49283;

dissociationEnergy = -(eLCAO[rMin a0] - eLCAO[20 a0]);

dissociationEnergy/e

eLCAOrescaled[r_] := eLCAO[r] + dissociationEnergy - e1

vMorse[dEq_, a_, rEq_, r_] := dEq (1 - Exp[-a (r - rEq)])^2

Plot[{vMorse[dissociationEnergy, 0.5, rMin, r]/e,vMorse[dissociationEnergy, 1.0, rMin, r]/e,eLCAOrescaled[r a0]/e}, {r,0, 12},AxesLabel->{"R/(Subscript[a, 0])","Energy(rescaled)/eV"}]

Plot[{vMorse[dissociationEnergy, 0.73, rMin, r]/e,eLCAOrescaled[r a0]/e}, {r,0, 12},AxesLabel->{"R/(Subscript[a, 0])","Energy(rescaled)/eV"}]