We present the extension of the complete flux scheme to advection-diffusion-reaction systems. The flux approximation is derived from a local system boundary value problem for the entire system, including the source term vector. Therefore, the numerical flux vector consists of a homogeneous and an inhomogeneous component, corresponding to the homogeneous and the particular solution of the boundary value problem, respectively. The complete flux scheme is validated for a test problem and shows uniform second order convergence behaviour