Dynamic thermography has been clinically proven to be a valuable diagnostic technique for skin tumor detection as well as for other medical applications, and shows many advantages over static thermography. Numerical modelling of heat transfer phenomena in biological tissue during dynamic thermography can aid the technique by improving process parameters or by estimating unknown tissue parameters based on measurement data. This paper presents a new non-linear numerical model of multilayer skin tissue containing a skin tumour together with thermoregulation response of the tissue during the cooling-rewarming process of dynamic thermography. The thermoregulation response is modelled by temperature-dependent blood perfusion rate and metabolic heat generation. The aim is to describe bioheat transfer more realistically. The model is based on the Pennes bioheat equation and solved numerically using a subdomain BEM approach treating the problem as axisymmetrical. The paper includes computational tests for Clark II and Clark IV tumours, comparing the models using constant and temperature-dependent properties which showed noticeable differences and highlighted the importance of using a local thermoregulation model. Results also show the advantage of using dynamic thermography for skin tumour screening and detection at an early stage. One of the contributions of this paper is a complete sensitivity analysis of 56 model parameters based on the gradient of the surface temperature difference between tumour and healthy skin. The analysis shows that size of the tumour, blood perfusion rate, thermoregulation coefficient of the tumour, body core temperature and density and specific heat of the skin layers in which the tumour is embedded are important for modelling the problem, and so have to be determined more accurately to reflect realistic skin response of the investigated tissue, while metabolic heat generation and its thermoregulation are not.