We propose in this paper an anisotropic, adaptive, finite element algorithm for steady, linear advection-diffusion-reaction problems with strong anisotropic features. The error analysis is based on the dual weighted residual methodology, allowing us to perform “goal-oriented” adaptation of a certain functional of the solution and derive an “optimal” metric tensor for local mesh adaptation with linear and quadratic finite elements. As a novelty, and to evaluate the weights of the error estimator on unstructured meshes composed of anisotropic triangles, we make use of a patchwise, higher-order interpolation recovery readily extendable to finite elements of arbitrary order. We carry out a number of numerical experiments in two dimensions so as to prove the capabilities of the goal-oriented adaptive method. We compute the convergence rate and the effectivity index for a series of output functionals of the solution. The results show the good performance of the algorithm with linear as well as quadratic finite elements.