% to plot 1sa orbital for H2+ for zeta=1
syms x y
% internuclear distance
R=4;
zeta=1;
% ra is the distance of the electron from the two nuclei
ra=sqrt(y^2+(x+R/2)^2);
% 1sa  orbital
onesa=sqrt(zeta^3/pi)*exp(-zeta*ra);

% contour plot of the density
ezcontour(onesa.*onesa,[-4 4 -4 4])
axis('square');
set(gca,'fontsize',14);title(['1s_a^2; R = ' num2str(R) '; zeta = ' num2str(zeta)])
set(gca,'xtick',[-4:4],'ytick',[-4:4])
hold;plot([-R/2 R/2],[0 0],'-k','linewidth',2);hold
