Detailed numerical calculations are used to evaluate the orientational distribution function of monodomain polymeric liquid crystals in simple shear and elongational flow. Doi's theory of concentrated solutions of rigid rods is used to model monodomain systems of nematic PLCs. In addition, an analysis of the validity of the closure approximation in Doi's theory is presented. A novel numerical technique is presented for the solution of the distribution function. This technique can also be applied to the evaluation of other distribution functions encountered in polymer or fiber suspension rheology. Results from the numerical calculations are presented as a function of the rod concentration for different flow kinematics. Qualitatively, agreement between the decoupled elastic stress and the exact values is good, but quantitatively, there is some discrepancy.