OK, having tried it again using a single Vcross term, I get the energy eigenvalues as:
E = 1/2 * [E1 + E2 + 2V +- sqrt ( (E1 - E2)^2 + 4|t|^2 ) ]
Where E1 is the epsilon01 for nucleus 1, sim. for E2, V is Vcross, t is the hopping term.
This makes sense as an expression, as the bit outside the square root is a sort of average of the two energies, and the difference between the energy eigenvalues depends on the difference between E1 and E2.
I would guess that to show the bonding orbital becomes more localised on the lower energy atom, it would be required to find the eigenvector for the lower eigenvalue... which I have done. The result isn't pretty, but it depends only on the difference E1-E2, which I would guess is a good thing.
If the V term is dropped (as in the "standard case" by a Gauss law type argument), then I get:
t*phi2 = [(E1-E2)/2 - sqrt(((E1-E2)/2)^2 + |t|^2)]*phi1
Here phi1 and phi2 have the same meaning as they do in the "standard case" presented in the main text.
Should it be obvious from this that the wavefunction becomes more localised on phi1 if E1 is lower (and vice versa)? (Also, does that expression look reasonable?)