How does one find the probability of tunneling along a reaction coordinate computationally?
If the reaction coordinate $\xi$ can be expressed as the change in the relative position between two groups of the molecule (e.g. for instance as a change in a bond length so that $\xi$ is given in units of length) the problem may be treated as one dimensional and one may use computational software to calculate the potential energy $V(\xi)$ along this path. You may then proceed solving the time-independent Schroedinger equation using numerical techniques such as Numerov's integration method. However, in many cases WKB theory provides an accurate estimate. 

In terms of WKB, the tunneling probability is given by
$$
T=e^{-2G}
$$
where $G$ is the Gamow factor given by
$$
G=\frac{2}{\hbar}\int_{\xi_a}^{\xi_b}{ \sqrt{2\mu_\text{red}\left [V(\xi) - E \right ]}}d\xi
$$
where $\mu_\text{red}$ is the reduced mass along the reaction coordinate and $E$ is the energy of the incoming particles. The integral runs over the clasically forbidden region from $\xi_a$ to $\xi_b$. The Gamow factor thus represents the square root of the integral over the "reaction barrier" multiplied with the effective mass that tunnels. Of course, you should take zero point energy into account and in the case of light reactants also the vibrational and rotational structure of the products as this may lead to resonances.