(* Lennard-Jones potential with a centrifugal term *)
a= 1.0; (* Depth of LJ well *)
l =1.0; (* Angular momentum *)
v12[r_] = 4 * a /r^12;  (* Repulsive hard-core potential *)
v6[r_] = -4 * a / r^6; (* Van der Waals potential *)
vCentr[r_]=0.5* l^2/r^2; (* Centrifugal potential *)
t[p_] = 0.5 * p^2; (* Kinetic energy *)
vLJ[r_]=v12[r]+v6[r]; (* Lennard-Jones potential *)
v[r_]=vLJ[r]+vCentr[r]; (* Effective potential *)
f[r_]=-D[v[r],r]; (* Force *)
h[r_]=-D[f[r],r]; (* Hessian *)
rMin = r/.NSolve[{f[r]==0, h[r]>0, r>0},r, Reals][[1]]
rMax = r/.NSolve[{f[r]==0, h[r]<0, r>0},r, Reals][[1]]
vMin = v[rMin];
vMax=v[rMax];

(* Purely repulsive r12 potential *)
Plot[v12[r], {r,0,3}, PlotRange->{-a,2}]
h[r_,p_]=t[p]+v12[r];
ContourPlot[h[r,p]=={0.5,1.,2., 3., 4.},{r,0,3},{p,-3,3}]

(* Lennard-Jones potential *) 
Plot[vLJ[r], {r,0,3}, PlotRange->{-a,2}]
h[r_,p_]=t[p]+vLJ[r];
ContourPlot[h[r,p]=={-0.98,-0.5 ,0.0, 0.5 , 1. ,2.}*a,{r,0,3},{p,-3,3}]

(* Lennard-Jones with centrifugal repulsion *)
h[r_,p_]=t[p]+v[r];
Plot[v[r], {r,0,1.55*rMax}, PlotRange->{-a,2}]
ContourPlot[h[r,p]=={0.98*vMin,0.5*vMin,0.0, 0.8*vMax,vMax*1.01,0.4,1.,2.},{r,0,1.5*rMax},{p,-3,3}]