Abstract We derived a continuous Galerkin spectral element method to compute approximations of the anisotropic diffusion operator with the tensorial coefficient of order two. After writing the operator in the weak form, we calculated spatial integration using the Legendre base. Then, we applied three different methods for the solution of the associated linear system: The Lower-Upper factorization, the biconjugate gradient method and the biconjugate gradient stabilized method with the Richardson preconditioner. To validate the algorithm, we present a convergence study for the Poisson equation when varying the functions of the tensor coefficient. Results show how the algorithm can be used to construct solvers for more complex problems like the anisotropic diffusion-reaction equation.