Diffusion of impurities in a granular gas undergoing homogeneous cooling state is studied. The results are obtained by solving the Boltzmann-Lorentz equation by means of the Chapman-Enskog method. In the first order in the density gradient of impurities, the diffusion coefficient D is determined as the solution of a linear integral equation which is approximately solved by making an expansion in Sonine polynomials. In this paper, we evaluate D up to the second order in the Sonine expansion and get explicit expressions for D in terms of the coefficients of restitution for the impurity-gas and gas-gas collisions as well as the ratios of mass and particle sizes. To check the reliability of the Sonine polynomial solution, analytical results are compared with those obtained from numerical solutions of the Boltzmann equation by means of the direct simulation Monte Carlo method. In the simulations, the diffusion coefficient is measured via the mean-square displacement of impurities. The comparison between theory and simulation shows in general an excellent agreement, except for the cases in which the gas particles are much heavier and/or much larger than impurities. In these cases, the second Sonine approximation to D improves significantly the qualitative predictions made from the first Sonine approximation. A discussion on the convergence of the Sonine polynomial expansion is also carried out.