We propose and analyze the ReLPM (Real Leja Points Method) for evaluating the propagator φ(ΔtB)v via matrix interpolation polynomials at spectral Leja sequences. Here B is the large, sparse, nonsymmetric matrix arising from stable 2D or 3D finite-difference discretization of linear advection-diffusion equations, and φ(z) is the entire function φ(z)=(e z -1)/z. The corresponding stiff differential system y(t)=By(t)+g,y(0)=y 0 , is solved by the exact time marching scheme y i + 1 =y i +Δt i φ(Δt i B)(By i +g), i=0,1,..., where the time-step is controlled simply via the variation percentage of the solution, and can be large. Numerical tests show substantial speed-ups (up to one order of magnitude) with respect to a classical variable step-size Crank-Nicolson solver.