The Legendre Galerkin Chebyshev collocation least squares method is presented for a second‐order elliptic problem with variable coefficients. By introducing a flux variable, the original problem is rewritten as an equivalent first‐order system. The present method is based on the Legendre Galerkin method, but Chebyshev–Gauss–Lobatto collocation is used to deal with the variable coefficient and the right hand side terms. The coercivity and continuity of the method are proved and an error estimate in the
H
1
‐norm is derived. Some numerical examples are given to validate the efficiency and accuracy of the scheme. © 2016Wiley Periodicals, Inc. Numer Methods Partial Differential Eq 32: 1689–1703, 2016