In this paper, we develop a density evolution (DE) framework for analyzing the iterative joint decoding (JD) for non-orthogonal multiple access (NOMA) systems, where the ordered-statistics decoding (OSD) is applied to decode short block codes. We first investigate the density-transform feature of the soft-output OSD (SOSD), by deriving the density of the extrinsic log-likelihood ratio (LLR) with known densities of the priori LLR. Then, we represent the OSD-based JD by bipartite graphs (BGs), and develop the DE framework by characterizing the density-transform features of nodes over the BG. Numerical examples show that the proposed DE framework accurately tracks the evolution of LLRs during the iterative decoding, especially at moderate-to-high SNRs. Based on the DE framework, we further analyze the BER performance of the OSD-based JD, and the convergence points of the two-user and equal-power systems.