Anderson acceleration method is an algorithm for accelerating the convergence of general fixed point iteration or Picard iteration. The method was first proposed in 1965 to accelerate the convergence of self-consistent field iteration in electronic structure computation, and it has being used successfully in that area ever since. However, until recently, Anderson acceleration method began to attract attention from the community of numerical analysis and algorithm development and some other application areas. Compared with Newton-like method, one advantage of Anderson acceleration method is that there is no need to form the Jacobian matrix. Therefore, the method is easy to be implemented. In this paper, Anderson accelerated Picard method is employed to solve a kind strong nonlinear radiation diffusion equation, three temperature energy equations. To improve the robustness of Anderson acceleration method, two strategies are incorporated into Anderson acceleration method. One strategy is used so that the iteration satisfying physical constraint. Another strategy is about monitoring matrix condition number for the least squares problem in the implementation of Anderson acceleration so that numerical stability is guaranteed. Numerical results show that Anderson accelerated Picard method can solve three temperature energy equation effectively. Compared to no accelerated Picard method, at least half iteration numbers can be saved by using Anderson acceleration method. Comparison for Jacobian-free Newton-Krylov method to Picard method and Anderson acceleration method is given.