Left-invariant PDE-evolutions on the roto-translation group $SE(2)$ (and their resolvent equations) have been widely studied in the fields of cortical modeling and image analysis. They include hypo-elliptic diffusion (for contour enhancement) proposed by Citti & Sarti, and Petitot, and they include the direction process (for contour completion) proposed by Mumford. This paper presents a thorough study and comparison of the many numerical approaches, which, remarkably, is missing in the literature. Existing numerical approaches can be classified into 3 categories: Finite difference methods, Fourier based methods (equivalent to $SE(2)$-Fourier methods), and stochastic methods (Monte Carlo simulations). There are also 3 types of exact solutions to the PDE-evolutions that were derived explicitly (in the spatial Fourier domain) in previous works by Duits and van Almsick in 2005. Here we provide an overview of these 3 types of exact solutions and explain how they relate to each of the 3 numerical approaches. We compute relative errors of all numerical approaches to the exact solutions, and the Fourier based methods show us the best performance with smallest relative errors. We also provide an improvement of Mathematica algorithms for evaluating Mathieu-functions, crucial in implementations of the exact solutions. Furthermore, we include an asymptotical analysis of the singularities within the kernels and we propose a probabilistic extension of underlying stochastic processes that overcomes the singular behavior in the origin of time-integrated kernels. Finally, we show retinal imaging applications of combining left-invariant PDE-evolutions with invertible orientation scores.