Coupling is a way of constructing Markov processes with prescribed laws on the same probability space. It is known that the rate of coupling (how fast you can make two processes meet) of elliptic/Riemannian diffusions is connected to the geometry of the underlying space. In this talk we consider coupling of hypoelliptic diffusions (diffusions driven by vector fields satisfying Hoermander's condition). S. Banerjee and W. Kendall constructed successful Markovian couplings for a large class of hypoelliptic diffusions. We use a non-Markovian coupling of Brownian motions on the Heisenberg group, and then use this coupling to prove analytic gradient estimates for harmonic functions for the sub-Laplacian. This talk is based on the joint work with Sayan Banerjee and Phanuel Mariano.