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 under the binary phase shift keying (BPSK) transmission. 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.