r = 0.5:0.05:5;
E = zeros(1,length(r));
for k = 1:length(r)
    E(k) = hf_potential(r(k)) + 2/r(k);
end

plot(r,E);
hold on;
plot(r,E,'o');
plot(r([5:30]),E([5:30]),'o');
grid on;
xlabel('He-H interatomic distance');
ylabel('The effective potential for nuclear motion');
title('The HF/Born-Oppenheimer approximation for effective potential of nuclear motion of HeH+');