Abstract: The spherical harmonics （P_N） method is applied to a second order Self-Adjoint Angular Flux （SAAF） neutron transport equation in the two-dimensional cylindrical coordinates. The spherical harmonics functions are used to expand the angular flux. A set of differential equations about the spatial variables is derived. The spatial variables are discretized by means of the finite element method. A computational code for multi-group neutron transport equation named TEPFEM-C is developed based on the present model. The numerical results of some benchmark problems demonstrate that this method can give high precision results and avoid the ray effect very well.