Polymorphism, in which a chemical compound exhibits different crystal forms or structures, has significant influence on the processing and storage of some crystalline powders in pharmaceutical industry. Different crystal structures, the so-called polymorphs, have different physical and chemical properties, such as crystal morphology, solubility, and color. These properties can have profound effect on the performance of products. This fact has motivated several researchers in this field to analyze, simulate, and control the crystallization of polymorphs. In this article, an efficient and accurate numerical method is introduced for solving a model describing crystallization of polymorphs. The proposed method has two parts. In the first part, a coupled system of ordinary differential equations of moments and the solute concentration is numerically solved in the time domain of interest. The resulting values are used to get the discrete growth and nucleation rates in the same time domain. In the second part, these discrete values along with the initial crystal size distribution (CSD) are used to construct the final CSD. The method is applied to a model describing crystallization of polymorphs of L-glutamic acid. For a validation, the results of the proposed technique are compared with those from the finite volume schemes. The numerical results demonstrate the potential of our scheme for the simulation of current model with high efficiency and accuracy.