To my knowledge this hasn't been done in theory (although see Harriet Moser's thesis http://www.math.columbia.edu/~moser/). But it certainly has been done in practice by Jeff Weeks' program SnapPea. Note that $\mbox{Isom}(\mathbb{H}^3) = \mbox{PSL}(2, \mathbb{C}) = \Gamma$. So your source group $G = \pi_1(M^3)$ already has a very nice matrix group as a target. SnapPea assumes that the three-manifold $M$ is given as a triangulation. (I am sure that it is not easy to go backwards from a presentation of $G$ to a triangulation of $M$.) After tidying the triangulation (and drilling out a curve if necessary - but lets suppose that $M$ has a single torus boundary component) SnapPea gives all of the tetrahedra in the triangulation the same shape, namely that of the regular ideal tetrahedron. "Developing" in $\mathbb{H}^3$ turns shapes of tetrahedra into matrices, one per generator. Naturally this is not yet a representation of $G$ into $\Gamma$. The failure to be a representation is measured by the failure of the shapes to satisfy the Thurston gluing equations. SnapPea uses a multivariate version of Newton's method to find new shapes for the tetrahedra, hopefully converging to the discrete and faithful representation of $G$ into $\Gamma$.
Details, references, and more can be found in Moser's thesis. However, I'll add a final remark - the convergence properties of SnapPea's method certainly do depend sensitively on the initial triangulation. There are certain triangulations where SnapPea will consistently produce wrong answers. This is not a problem in practice -- you randomize the triangulation a few times and SnapPea typically starts to behave much better. But finding an actual algorithm appears to be difficult. So something mysterious going on.
The problem (of finding the discrete and faithful representation) is easy to solve even for very respectfully sized manifolds (say up to 100 tetrahedra). But why? I certainly don't know. You'd have to ask Thurston, Weeks, Hodgson or some other expect expert for an opinion.
