Methods for determining and computing the rate-distortion (RD) bound for N-layer scalable source coding of a finite memoryless source are considered. Optimality conditions were previously derived for two layers in terms of the reproduction distributions q(y1) and q(y2/y1). However, the ignored and seemingly insignificant boundary cases, where q(y1) = 0 and q(y2/y1) is undefined, have major implications on the solution and its practical application. We demonstrate that, once the gap is filled and the result is extended to N-layers, it is, in general, impractical to validate a tentative solution, as one has to verify the conditions for all conceivable q(yi+1,...,yN/y1,...,yi) at each (y(1), ..., y(i)) such that q(y1,...,yi) = 0. As an alternative computational approach, we propose an iterative algorithm that converges to the optimal joint reproduction distribution q(y1',...,yN), if initialized with q(y1,...,yN) > 0 everywhere. For nonscalable coding (N = 1), the algorithm specializes to the Blahut-Arimoto algorithm. The algorithm may be used to directly compute the RD bound, or as an optimality testing procedure by applying it to a perturbed tentative solution q. We address two additional difficulties due to the higher dimensionality of the RD surface in the scalable (N > 1) case; namely, identifying the sufficient set of Lagrangian parameters to span the entire RD bound; and the problem of efficient navigation on the RD surface to compute a particular RD point.