In this paper, we propose a mixed Fourier-Jacobi spectral method for two dimensional Neumann boundary value problem. This method differs from the classical spectral method. The homogeneous Neumann boundary condition is satisfied exactly. Moreover, a tridiagonal matrix is employed, instead of the full stiffness matrix encountered in the classical variational formulation. For analyzing the numerical error, we establish the mixed Fourier-Jacobi orthogonal approximation. The convergence of proposed scheme is proved. Numerical results demonstrate the efficiency of this approach.