ABSTRACT:The pairs of radial functions P i and Q i , which are part of the four-component single-particle spinors in the relativistic description of the electronic structure of bound states of atoms, are usually determined as solutions of eigenvalue problems. The latter constitute two-point boundary value problems which involve coupled pairs of first-order ordinary radial differential equations. To introduce a suitable notation, the theory involved in relativistic electronic structure calculations for atoms is briefly reviewed, including a general treatment of arbitrary transformations for the radial variable. Such variable transformations must be specified then, either explicitly (by a transformation function) or implicitly (by the solution method itself), for any actual calculation, though initially the variable transformation is almost completely independent of the numerical method to be applied. We consider suitable transformation functions out of which an actual choice can be made within wide limits. In our approach, the resulting transformed radial differential equations are discretized then with finite-difference methods. Since the accuracy and efficiency of these methods is increased when a constant step size h between contiguous points (in the transformed variable) is used, it is only here where the careful choice of the variable transformation is important. The resulting system of linear equations is solved for the radial functions with standard linear algebra methods rather than "shooting" methods. The present work extends the general study of variable transformations, given recently for nonrelativistic electronic structure calculations in Part I, to the relativistic case. Important results following from the present work are (i) a discretization scheme for first-order ordinary differential equations similar to the ANDRAE, REIHER, AND HINZE well-known standard Numerov scheme used within the nonrelativistic framework, (ii) a consistent numerical algorithm with an overall truncation error of order h 4 , (iii) a method for handling effective potentials behaving singularly like r −1 at the origin (as encountered, e.g., when the atomic nucleus is represented by a pointlike charge density distribution), and, in connection with this last point, (iv) the avoidance of nonanalytic short-range behavior of the solutions to be obtained from the differential equations.