It turns out we already have everything we need to compute the turning points - and thus to compute the potential. To this end, the program loops through values of v from vmin (usually -0.5 unless we apply the Kaiser correction) to vmax in increments of space. This is fully accomplished by the above code snippet. As the turning points are calculated, they (and their associated energy) are added to an ever-expanding list.

The above figure presents the raw output of iteratively applying the RKR method on carbon monoxide up to a maximum vibrational level of 80 in increments of 2. Each cross represents a single computed turning point, while the lines connecting them are the result of simple linear interpolation. In any application of this program, we would use a much smaller value for 'space' (so the computed turning points are much denser, which is more useful for subsequent calculations). However, it is instructive to use 2 here (so the discrete turning points can be directly observed).

There are a couple of problems here: on the left-hand side, where the potential should go to infinity (see the molecular potential diagram on the homepage), it instead doubles back - a known problem with the RKR method, as the approximations used to make it valid break down for high vibrational levels. A little trickier to spot is that the dissociation energy is too high - the potential should level out much sooner than we observe here. So, some corrections are in order.

The program performs a spline interpolation on the (ideally very dense) array of turning points, producing a highly accurate piecewise function which "fills in" the points on the potential curve that we didn't explicitly calculate.

Using the resultant spline object, as well as its first and second derivatives, the program performs a simple search to locate downwards curvature on the inner portion of the curve and the location of the potential minimum.

If the program is able to identify negative concavity before the well minimum,  it recognizes that the raw RKR data is in need of a correction. Previous work by Robert J. LeRoy suggested that an inverse-power fit of the form (A/r^n + B) suffices for this kind of correction. To obtain this fit, the program computes the best value of n for each turning point near the problematic region, averages them out, then selects A and B to overlap the last two valid turning points.

The first equation above gives the best value for n for a given r. The second and third give the exact values for A and B given that v1 and v2 are the last two valid vibrational levels to be fitted to.