This paper deals with the problem of laminar axisymmetric incompressible flow in the working area of a rotary drum furnace. A system of second-order partial differential equations is derived and solved numerically using the finite difference method. The obtained solutions are tested for mesh convergence and compared with available analytical solutions for a laminar axisymmetric flow. The numerical results show that furnace rotation significantly affects the gas flow in the working section of the apparatus. Thus, the temperature distribution over the working area and furnace surface can be controlled through the furnace rotation velocity and gas flow rate. The input velocity profile also affects the flow pattern in the furnace. Vortex zones with counterflow occur near the walls, where the axial velocity profile is more elongated than the Poiseuille parabola. On the symmetry axis, the axial velocity is higher than that of the Poiseuille flow. The proposed mathematical model allows one to calculate more accurately the degree of thermal decomposition as a function of temperature and longitudinal coordinate along the furnace axis.