Modification of thin metallic films using ultrashort laser pulses involves interplay of numerous physical processes. Finding a right combination of laser parameters is essential for achieving the desired modification of a thin film deposited on a substrate. Numerical modeling is a convenient tool for gaining insights into ultrafast evolution of material properties and to predict an optimal range of irradiation parameters. In this work, a mathematical model is presented that describes the ultrafast laser heating and temperature relaxation in a thin molybdenum film deposited on a glass substrate. The laser energy absorption by molybdenum is described using a two-temperature model. The model takes into account the heat exchange between the film and the substrate through a boundary condition applied on the lattice temperature. The implicit numerical scheme employed for simulations was verified in respect of energy conservation. The model has been validated by comparison with experimental data on melting threshold fluences.