rr=[.001 .1:.1:10];
global R
for i=1:length(rr);
   R=rr(i);
   zet(i)=fminsearch('h2plus_varzeta',1);
   e(i)=h2plus_varzeta(zet(i));
   ee(i)=h2plus_varzeta(1);
end
figure(1)
plot(rr,e,rr,ee)
set(gca,'fontsize',14);xlabel('R / bohr');ylabel('Electronic energy / hartree');
title('H_2^+; electronic energy; blue:  variable \zeta; green:  \zeta=1')
axis([0.1 10 -2 -.5])

figure(2)
plot(rr,zet)
set(gca,'fontsize',14);xlabel('R / bohr');ylabel('\zeta');
hold;plot([min(rr) max(rr)],[1 1],'-k');hold
title('H_2^+; electronic energy; best value of \zeta')


%%
[ee mn]=min(e+1./rr)
cc=polyfit(rr(mn-2:mn+2),e(mn-2:mn+2)+1./rr(mn-2:mn+2),4);
roots(polyder(cc));
R=ans(3);
fminsearch('h2plus_varzeta',1);
emin=h2plus_varzeta(ans)+1./R;
rmin=R;
de=-0.5 -(emin+1/R);
figure(3)
plot(rr,e+1./rr)
axis([.5 6 -0.605 -.5])
set(gca,'fontsize',14);xlabel('R / bohr');ylabel('V(R) / hartree')


% exact values from V. A. Johnson, Phys. Rev. 60, 373 (1941)

R=[1.214301 1.366077 1.523336 1.685909 1.853595 1.997466 2.026162 2.203349 2.384868 2.570407 2.759633
  ];
E=-[.065698 .129506 .169259 .192061 .202863 .205268 .205184 .201567 .193872 .183471 .171381]-1;
% convert from rydbergs to hartree
E=E/2;
hold;plot(R,E,'o');hold
title('H_2^+ Potential Curve; circles are exact values')

