A script that works for arbitrary potentials
My script has produced the following figures for this paper.
Simple single-well potential:

Double-minimum potential with a small barrier and a "shelf":

Notice that when the potential gets wider, the vibrational levels get closer together, as we expect based on what happens to a "particle in a box" when the box gets wider.
The same basic script was used for this paper.
Distinction among vibrational levels that have been spectroscopically observed versus ones that have not been observed:

Inset showing vibrational levels:

Vibrational levels labeled:
The same basic script was also used for the diagram in my MMSE answer here.

It was also used for the figure in this paper.
Selected vibrational levels are shown or not shown:

Here is the script that was used for the first paper listed above
close('all')
figure1=figure(1);
set(gca,'Position',[0.101190476190476 0.105359342915811 0.885 0.888357289527721])
set(gcf,'Color','w')
minX=1;maxX=20;
splineMeshSize_for_r=0.00001; % I just chose it to be of the same precision as r_e defined in cell 1. Splines can only be necessary for short-range, since long-range seems to have points very close together.
splineMesh_for_r=minX:splineMeshSize_for_r:maxX;
V_ab_initio_spline_6X=spline(r_abInitio_6X,V_abInitio_6X,splineMesh_for_r);
V_ab_initio_spline_6a=spline(r_abInitio_6a,V_abInitio_6a,splineMesh_for_r);
V_ab_initio_spline_3A=spline(r_abInitio_3A,V_abInitio_3A,splineMesh_for_r);
V_ab_initio_spline_3B=spline(r_abInitio_3B,V_abInitio_3B,splineMesh_for_r);
V_ab_initio_spline_3C=spline(r_abInitio_3C,V_abInitio_3C,splineMesh_for_r);
V_ab_initio_spline_3b=spline(r_abInitio_3b,V_abInitio_3b,splineMesh_for_r);
V_ab_initio_spline_3c=spline(r_abInitio_3c,V_abInitio_3c,splineMesh_for_r);
V_ab_initio_spline_3d=spline(r_abInitio_3d,V_abInitio_3d,splineMesh_for_r);
minY=-6000;maxY=300;
axis([minX,maxX,minY,maxY])
line([minX maxX], [0 0],'Color','k','LineWidth',3)
hold('on')
% [Vmin_ab_initio, indexOfVmin_ab_initio]=min(V_ab_initio_spline);
% plotPointsSymmetricallyOnPotential(vibrationalEnergies(1:2:end),V_ab_initio_spline,splineMesh_for_r,indexOfVmin_ab_initio);
axis([minX,maxX,minY,maxY])
plot(r_abInitio_6X,V_abInitio_6X,'Color',[0,255,100]./255,'Linewidth',3);
plot(r_abInitio_3A,V_abInitio_3A,'Color','r','Linewidth',3);
plot(r_abInitio_3B,V_abInitio_3B,'Color','k','Linewidth',3);
plot(r_abInitio_3C,V_abInitio_3C,'Color','c','Linewidth',3);
figure2=figure(2);
set(gca,'Position',[0.101190476190476 0.105359342915811 0.885 0.888357289527721])
set(gcf,'Color','w')
minX=1;maxX=20;
minY=-10000;maxY=300;
axis([minX,maxX,minY,maxY])
line([minX maxX], [0 0],'Color','k','LineWidth',3)
splineMeshSize_for_r=0.00001; % I just chose it to be of the same precision as r_e defined in cell 1. Splines can only be necessary for short-range, since long-range seems to have points very close together.
splineMesh_for_r=minX:splineMeshSize_for_r:maxX;
hold('on')
plot(r_abInitio_6a,V_abInitio_6a,'Color','b','Linewidth',3);
plot(r_abInitio_3b,V_abInitio_3b,'Color','m','Linewidth',3);
plot(r_abInitio_3c,V_abInitio_3c,'Color','g','Linewidth',3);
plot(r_abInitio_3d,V_abInitio_3d,'Color',[100,100,100]./255,'Linewidth',3);
% [Vmin_ab_initio, indexOfVmin_ab_initio]=min(V_ab_initio_spline);
% plotPointsSymmetricallyOnPotential(vibrationalEnergies(1:2:end),V_ab_initio_spline,splineMesh_for_r,indexOfVmin_ab_initio);
axis([minX,maxX,minY,maxY])
plot(r_abInitio_6a,V_abInitio_6a,'Color','b','Linewidth',3);
plot(r_abInitio_3b,V_abInitio_3b,'Color','m','Linewidth',3);
plot(r_abInitio_3c,V_abInitio_3c,'Color','g','Linewidth',3);
plot(r_abInitio_3d,V_abInitio_3d,'Color',[100,100,100]./255,'Linewidth',3);
%%
close('all');figure3=figure(3);hold('on')
set(gca,'Position',[0.101190476190476 0.114161849710983 0.892953008436171 0.879554782732549])
set(gcf,'Color','w')
minX=2;maxX=15;
minY=-6000;maxY=300;
axis([minX,maxX,minY,maxY])
line([minX maxX], [0 0],'Color','k','LineWidth',3)
re=3.98;
plot(r_abInitio_3b,V_abInitio_3b,'Color',[0,255,100]./255,'Linewidth',10);
plot(r_3b_p6q7r05_90,V_3b_p6q7r05_90,'k','LineWidth',3)
v0=-5654.932449854602055;
v( 1)=-5493.765135128664951;
v( 2)=-5385.710984514686970;
v( 3)=-5241.269953247034209;
v( 4)=-5083.976206177168933;
v( 5)=-4914.85;
v( 6)=-4737.691247106215087;
v( 7)=-4554.829936403927604;
v( 8)=-4368.240502184650722;
v( 9)=-4179.292096277415112;
v(10)=-3988.756231213321371;
v(11)=-3796.972235927646125;
v(12)=-3604.218750080568952;
v(13)=-3411.033149015913750;
v(14)=-3218.297919388399805;
v(15)=-3027.108781656958399;
v(16)=-2838.477232599816034;
v(17)=-2652.930601841349471;
v(18)=-2470.302596380747673;
v(19)=-2290.105676375571875;
v(20)=-2112.256988503596288;
v(21)=-1937.405398146945536;
v(22)=-1766.631273238948864;
v(23)=-1600.818435795562496;
v(24)=-1440.124051819757312;
v(25)=-1284.340696145073408;
v(26)=-1134.583144744511360;
v(27)= -996.778888028432000;
v(28)= -914.777113202062848;
v(29)= -882.576799450836608;
v(30)= -845.026425674076544;
v(31)= -800.327833982061696;
v(32)= -753.924510206392448;
v(33)= -705.900692148498688;
v(34)= -654.935744351130112;
v(35)= -602.077023609842304;
v(36)= -548.234951970257472;
v(37)= -493.706925401166272;
v(38)= -438.926074728071424;
v(39)= -384.474432243322304;
v(40)= -330.903787772794240;
v(41)= -278.814572462308608;
v(42)= -228.913912424278848;
v(43)= -182.018345069401152;
v(44)= -139.060176467093648;
v(45)= -101.066168262798576;
v(46)= -69.047731355521066;
v(47)= -43.749645705548336;
v(48)= -25.299916111188852;
v(49)= -13.022995406988458;
v(50)= -5.669747141689355;
v(51)= -1.855271158639904;
v(52)= -0.324750667836605;
v(53)= -0.002393731491646;
vibrationalEnergies=[v0 v];
[Vmin_ab_initio_3b, indexOfVmin_ab_initio_3b]=min(V_ab_initio_spline_3b);
%plotPointsSymmetricallyOnPotential(vibrationalEnergies(1:2:end),V_ab_initio_spline_3b,splineMesh_for_r,indexOfVmin_ab_initio_3b,200);
vibrationalLevelsMeasured=[0:53];
for vibrationalLevel=vibrationalLevelsMeasured+1
line([interp1(V_3b_p6q7r05_90(r_3b_p6q7r05_90<re),r_3b_p6q7r05_90(r_3b_p6q7r05_90<re),vibrationalEnergies(vibrationalLevel)) interp1(V_3b_p6q7r05_90(r_3b_p6q7r05_90>re),r_3b_p6q7r05_90(r_3b_p6q7r05_90>re),vibrationalEnergies(vibrationalLevel))], [vibrationalEnergies(vibrationalLevel) vibrationalEnergies(vibrationalLevel)],'Color','b','LineWidth',3)
r_innerTurningPoint(vibrationalLevel)=interp1(V_3b_p6q7r05_90(r_3b_p6q7r05_90<re),r_3b_p6q7r05_90(r_3b_p6q7r05_90<re),vibrationalEnergies(vibrationalLevel));
scatter(r_innerTurningPoint(vibrationalLevel),vibrationalEnergies(vibrationalLevel),'MarkerFaceColor','b','MarkerEdgeColor','b');
r_outerTurningPoint(vibrationalLevel)=interp1(V_3b_p6q7r05_90(r_3b_p6q7r05_90>re),r_3b_p6q7r05_90(r_3b_p6q7r05_90>re),vibrationalEnergies(vibrationalLevel));
scatter(r_outerTurningPoint(vibrationalLevel),vibrationalEnergies(vibrationalLevel),'MarkerFaceColor','b','MarkerEdgeColor','b');
end
plot(r(700:end),V(end,700:end),'r','LineWidth',6)
annotation(figure3,'textbox',[0.620 0.760 0.327 0.087],'Color','r','String','$V(r) = - \frac{D_3(r)C_3}{r^3} - \frac{D_6(r)C_6}{r^6} - \frac{D_8C_8(r)}{r^8}$','LineStyle','none','Interpreter','latex','FontSize',24);
annotation(figure3,'textbox',[0.421 0.609 0.107 0.0636],'String','MLR$_{5.9}^{6,7}(17)$','LineStyle','none','Interpreter','latex','FontSize',24);
annotation(figure3,'textbox',[0.293 0.223 0.107 0.0636],'String','Li$_2\left(3b,3^3\Pi_u\right)$','LineStyle','none','Interpreter','latex','FontSize',36,'FontName','Helvetica','FitBoxToText','off');
annotation(figure3,'textbox',[0.438 0.532 0.107 0.0636],'Color',[0,255,100]./255,'String','Original','LineStyle','none','Interpreter','latex','FontSize',24);
annotation(figure3,'line',[0.366764275256223 0.412884333821376],[0.63775968992248 0.637209302325581],'LineWidth',3);
scatter(5.89,-2750,200,'MarkerFaceColor',[0,255,100]./255,'MarkerEdgeColor','k')
scatter(6.5,-2750,200,'MarkerFaceColor',[0,255,100]./255,'MarkerEdgeColor','k')
annotation(figure3,'arrow',[0.620790629575403 0.588579795021962],[0.850162790697674 0.924031007751938],'LineWidth',3);
annotation(figure3,'arrow',[0.922401171303075 0.953147877013177],[0.768992248062015 0.689922480620155],'LineWidth',3);
xlabel('Internuclear distance \AA','Interpreter','Latex','FontSize',36) % (r) taken out as per request by Kirk Madison in email on 25/8/2012
ylabel('$V(r)$ cm$^{-1}$','Interpreter','Latex','FontSize',36)
box('on');
set(gca,'XMinorTick','on','YMinorTick','on','LineWidth',2,'FontSize',20);