The problem of the computation of the Centrifugal Distortion Constants (CDC) related to a diatomic potential is considered. The analytical expressions obtained from a reformulation of the Rayleigh-Schrodinger perturbation theory are used [Kobeissi et al., J. Mol. Spectrosc., 138, 1 (1989)]; these are e(n+1) = [PHI0RPHI(n)] -SIGMA(m=1)n e(m)[PHI0PHI(n-m)] where R = 1/r2, PHI0 = psi(v) is the vibrational wave function (corresponding to the given energy E(v) = e0) and PHI1, PHI2, . . . , are the "rotational corrections" to PHI0, solutions of the rotational (nonhomogeneous) Schrodinger equations. These equations are integrated by using a recent integrator using a powerful local control allowing (for PHI0) a high accuracy. The integrals are computed by using another powerful technique tailored for matrix elements between numerical wave functions [Kobeissi et al., J. Comp. Chem., 10, 358 (1989)]. This numerical treatment is applied to the model Lennard-Jones potential and to the RKR potential of the I2 ground state. In both applications the CDC are computed up to e6 = N(v) and e7 = O(v) (these two are published for the first time), and up to the dissociation [up to v = 23 for the Lennard-Jones potential, and to v = 108 for the XSIGMA - I2 (RKR) potential].