ABSTRACTIn this paper, we investigate the Legendre spectral methods for problems with the essential imposition of Neumann boundary condition in three dimensions. A double diagonalization process has been employed, instead of the full stiffness matrices encountered in the classical variational formulation of the problem with a weak natural imposition of Neumann boundary condition. For analysing numerical errors, some results on Legendre orthogonal approximation in Jacobi weighted Sobolev space are established. As examples of applications, the spectral schemes are provided for two model problems. The convergences of the proposed schemes are proved, too. Numerical results demonstrate the spectral accuracy in space, and which confirm theoretical analysis well.