Cebeci \& Bradshaw``s improved mixing length model, Mellor \& Herring``s turbulent kinetic energy equation model, and Ng \& Spalding``s two equation model are examined with wall bounded shear flows. Two dimensional flat plate boundary layer flows and fully developed channel flows are selected as test flows. An implicit finite difference scheme with variable grid size is used to solve the modeled equations. Since the two equation model is not valid for very low turbulent Reynolds number region, a computational matching scheme is developed in this study in such a way that within an equilibrium layer solutions of the two equation model are smoothly matched with those of Mellor \& Herring``s model. By comparing various models, it is found that the Ng \& Spalding``s two equation model with matching boundary conditions yields best results.