A mathematical model of the blast furnace was developed and the effects of operating variables on internal phenomena of the blast furnace were investigated through numerical computation. To achieve these objectives, equations which represent the distributions of gas and solid temperatures, gas compositions, pressure and stepwise reduction of iron ore were established and then solved simultaneously. To reduce computation time, orthogonal collocation in radial direction was employed after transforming the blast furnace geometry into a cylinder using the technique proposed by Yun et al.[1]. Besides these, a potential instability problem of the model was analyzed using a simplified 1-dimensional model. With the aid of this model, initial values of gas and solid temperature at the stock line, which give a meaningful stable solution, could be determined. Collocations with four, five and six internal points were tried for comparison but the results were similar. Simulation results based on the POSCO’s operating data showed good agreement with the internal phenomena of real blast furnaces as well as with the results by other rescearchers. One of the important consequences through this work is that the temperature distributions are sensitively affected by a small change in the ratio of solid to gas flowrates.