Although heat treatment is a process of great technological importance in order to obtain desired mechanical properties such as hardness, the process was required a tedious and expensive experimentation to specify the process parameters. Consequently, the availability of reliable and efficient numerical simulation program would enable easy specification of process parameters to achieve desired microstructure and mechanical properties without defects like crack and distortion. In present work, the developed numerical simulation program could predict distributions of microstructure and thermal stress in steels under different cooling conditions. The computer program is based on the finite difference method for temperature analysis and microstructural changes and the finite element method for thermal stress analysis. Multi-phase decomposition model was used for description of diffusional austenite decompositions in low alloy steels during cooling after austenitization. The model predicts the progress of ferrite, pearlite, and bainite transformations simultaneously during quenching and estimates the amount of martensite also by using Koistinen and Marburger equation. To verify the developed program, the calculated results are compared with experimental ones of casting product. Based on these results, newly designed heat treatment process is proposed and it was proved to be effective for industry.