This paper aims to develop a numerical approximation for the solution of the advection–diffusion equation with constant and variable coefficients. We propose a numerical solution for the equation associated with Robin’s mixed boundary conditions perturbed with a small parameter, ε. The approximation is based on a couple of methods: a spectral method of Galerkin type with a basis composed of Legendre polynomials and a Gauss quadrature of type Gauss–Lobatto applied for integral calculations with stability and convergence analysis.

In addition, a Crank–Nicolson scheme is used for the temporal solution as a finite difference method. Several numerical examples are discussed to show the efficiency of the proposed numerical method, especially when ε tends to zero so that we obtain the exact solution of the classical problem with homogeneous Dirichlet boundary conditions. The numerical convergence is well presented in different examples. Therefore, we build an efficient numerical method for different types of partial differential equations with different boundary conditions.