The characteristics of crystallization in a batch Couette flow crystallizer are modeled by population balance, mass and energy equations. The effective and powerful method of Legendre polynomials is employed to transform the partial differential equation of population balance into a series of ordinary differential equations of Legendre series expansion coefficients. Optimization theory is applied to find the minimum coefficient of variations, which accounts for the uniformity of crystal‐size distribution within the crystallizer. The results are obtained by the Runge‐Kutta‐Gill computational method and are explained by the physical significance of the growth and nucleation phenomena.