In this work we analyze Thermal Desorption Spectroscopy (TDS) experiments of gases considering simultaneously atomic diffusion in the bulk and recombinative reaction on the surface of a solid sample. We perform computer simulations in a three-dimensional cubic lattice, and compare the simulated spectra with numerical results from a simple reaction-diffusion model, and analytical results from models for the fast diffusion and fast reaction regimes. We analyze the effects on the spectra of each relevant parameter: diffusivity, reaction coefficient, heating rate, initial concentration of particles and sample size. We also develop a rate-limiting step identification procedure based on the analysis of peak position as a function of initial concentration. This last procedure applied to the simulations and to experimental TDS spectra produces very good results.