The ice accreted on the airfoil is one of the critical drivers that causes the degradation of aerodynamic performance as well as aircraft accidents. Hence, an efficient numerical code to predict the accreted ice shape is crucial for the successful design of de-icing and anti-icing devices. To this end, a numerical code has been developed for the prediction of glaze ice accretion shape on 2D airfoil. Constant Source-Doublet method is used for the purpose of computational efficiency and heat transfer in the icing process is accounted for by Messinger model. The computational results are thoroughly compared against available experiments and other computation codes such as LEWICE and TRAJICE. The direction and thickness of ice horn are shown to yield similar results compared to the experiments and other codes. In addition, the effects of various parameters - temperature, free-stream velocity, liquid water contents, and droplet diameter - on the ice shape are systematically analyzed through parametric studies.