Solving diffusion problems requires numerical methods able to capture the heterogeneity over complex geometries and are robust in terms of positivity preserving, nonlinearity, and radiation diffusion. Current deep learning methods, although mesh-free, encounter difficulties in achieving convergence and exhibit low accuracy when confronted with these specific issues. In this paper, we develop a novel method to overcome these issues based on the recently proposed Random Feature Method (RFM). Our contributions include: (1) for anisotropic and discontinuous coefficient problems, we rewrite a diffusion problem into a first-order system and construct the corresponding loss function and approximation spaces; (2) to avoid negative solutions, we employ the square function as the activation function to enforce the positivity and the trust-region least-square solver to solve the corresponding optimization problem; (3) for the radiation diffusion problem, we enrich the approximation space of random feature functions with the heat kernel. Numerous numerical experiments demonstrate that the current method outperforms the standard RFM as well as deep learning methods in terms of accuracy, efficiency, and positivity preserving.