%constants
E1 = -13.6;
R = 0.074;
a0 = 0.0529;

%matrix elements
R0 = R/a0;

a = 2*E1*(1-(1+R0)*exp(-2*R0))/R0;
b = 2*E1*(1+R0)*exp(-R0);
s = exp(-R0)*(1+R0+(R0^2/3));

%matrices
H_u = [E1 + a, E1*s+b; E1*s+b, E1 + a];
S_u = [1, s; s, 1];

%find eigenvalues and eigenvectors
[vectors,energies] = eig(inv(S_u)*H_u);
z = linspace(-2,2);

z_L = -0.37;
z_R = 0.37;
a0 = 0.529;

u_L = exp(-abs(z - z_L)./a0)./sqrt(pi*a0^3);

u_R = exp(-abs(z - z_R)./a0)./sqrt(pi*a0^3);

phi_B = u_L + u_R;

phi_A = u_L - u_R;

figure(1);

plot(z, phi_B.^2, 'k-');
hold on
plot(z, phi_A.^2, 'b-');
grid on;
xlabel('z [Angstroms]');
ylabel('probability [arbitrary scaling]');
legend('Bonding', 'Antibonding');