An algorithm is presented which solves the multi-dimensional advection-diffusion equation on complex shapes to second-order accuracy and is asymptotically stable in time. This bounded-error result is achieved by constructing, on a rectangular grid, a differentiation matrix whose symmetric part is negative definite. The differentiation matrix accounts for the Dirichlet boundary condition by imposing penalty-like terms. Numerical examples in two dimensions show that the method is effective even where standard schemes, stable by traditional definitions, fail. It gives accurate, non oscillatory results even when boundary layers are not resolved.