The RKR method determines the classical turning points $r_{\pm}$ of the potential energy based on knowing experimentally determined spectroscopic constants. ( Rydberg, 1931, Z. Physik. v73,376, Klein, 1932 Z. Physik. v76, 226, Rees Proc. Phys. Soc.(Lond.) 1947, v59, 998.)
The energy equation takes the usual form. Parameters all in cm$^{-1}$, distances in m, reduced mass $\mu$ in u, $h$ in cm$^{-1}$s.
$$E_v=\omega_e(v+1/2) - x_e\omega_e(v+1/2)^2 + y_e\omega_e(v+1/2)^3 + z_e\omega_e(v+1/2)^4 +\cdots$$
The vibrational quantum number $v$ should be considered as a variable that has values from $-1/2$ to whatever positive value is required.
The derivative of the rotational energy wrt. $J$ at $J=0$ is also needed,
$$B'(v)= B_e-\alpha_e(v+1/2)+\gamma_e(v+1/2)^2+\cdots$$
The function representing the turning points is
$$r_\pm= \frac{f(v)}{2}\left (\sqrt{1+\frac{4}{f(v)g(v)}}\pm 1 \right)$$
where functions $f, \, g$ are integrals. These functions 'blow up' at the upper integration limit $x \to v$ so a correction is added to prevent this. The basic equation is
$$f(v)=\frac{h}{\pi\sqrt{2\mu}}\int_{-1/2}^v\frac{dx}{ \sqrt{E(v)-E(x)}}$$
and the modified one
$$f(v)=\frac{h}{\pi\sqrt{2\mu}}\left(\int_{-1/2}^{v-\delta}\frac{dx}{ \sqrt{E(v)-E(x)} }+ 2\sqrt{\frac{\delta}{Q_v}}\right)$$
where $\delta$ is small $\approx 10^{-5}$ and
$$Q_v=\omega_e - 2x_e\omega_e(v+1/2) + 3y_e\omega_e(v+1/2)^2 + 4z_e\omega_e(v+1/2)^3 +\cdots$$
The function $g$ with correction is
$$g(v)=\frac{4\pi\sqrt{2\mu}}{h} \left( \int_{-1/2}^{v-\delta} \frac{B'(x)}{\sqrt{E(v)-E(x)} }dx+2B'(v)\sqrt{\frac{\delta}{Q_v} } \right) $$
A set of values for CO is
$\mu = 6.85841,\; \omega_e= 2169.8135,\; x_e\omega_e=13.2883,\; y_e\omega_e=0.010511,\; z_e\omega_e= 0.000057$
$ B_e= 1.9312808,\; \alpha_e= 0.0175044,\; \gamma_e= 0.000000548$.
Plot energy $E(x)$(cm$^{-1}) $ vs $r_\pm(x)$ (pm) is shown. The value at $x=-1/2$ should be $r_e=103.14$ pm the equilibrium bond length with energy of zero. As values are based on cm$^{-1}$ a scaling of $10^{14}$ is needed to go to pm.
I used Jupyter notebooks (via Anaconda) and python 3 to do the calculation using built in integrator from numpy/scipy and matplotlib to draw the graph. All are free and very easy to use.
