We present a hybrid method for the numerical solution of advection-diffusion problems that combines two standard algorithms: semi-Lagrangian schemes for hyperbolic advection-reaction problems and Crank-Nicolson schemes for purely diffusive problems. We show that the hybrid scheme is identical to the two end-member schemes in the limit of infinite and zero Peclet number and remains accurate over a wide range of Peclet numbers. This scheme does not have a CFL stability criterion allowing the choice of time step to be decoupled from the spatial resolution. We compare numerical results with an analytic solution and test both an operator split version of our method and a combined version that solves advection and diffusion simultaneously. We also compare results of simple explicit and implicit numerical schemes and show that the semi-Lagrangian Crank-Nicolson (SLCN) scheme is both faster and more accurate on the same problem. We then apply the combined SLCN scheme to a more geologically relevant benchmark for calculating the thermal structure of a subduction zone. This problem demonstrates that the SLCN scheme can remain stable and accurate at large Courant numbers even in flows with highly curved streamlines. Finally, we introduce a variable order interpolation scheme for the semi-Lagrangian schemes that reduces interpolation artifacts for sharp fronts without introducing numerical diffusion.