Yes! Most methods to compute exponentials of large sparse matrices are based on computing directly $\exp(A)b$ for a given vector $b$ rather than the full matrix $\exp(A)$. Just take $b$ as a vector of the canonical basis, and you're set.

The basic idea is that these algorithms are based on approximating the exponential with a rational function which is then expanded into a sum of [partial fractions](https://en.wikipedia.org/wiki/Partial_fraction_decomposition)
$$\exp(A)b \approx q(A)^{-1}p(A)b = \sum_{i=1}^k \omega_k (A-\tau_k I)^{-1} b.$$

Hence the problem is reduced to solving several shifted linear systems of the form $(A-\tau_k I)x_k=b$; this should work particularly well for a tridiagonal matrix, for which linear systems can be solved cheaply. 

This technique is often combined with ideas from [Krylov subspace](https://en.wikipedia.org/wiki/Krylov_subspace) methods, especially rational Krylov methods.

See for instance a chapter in the review  https://doi.org/10.1017/S0962492910000036,  https://doi.org/10.1137/100788860, or http://onlinelibrary.wiley.com/doi/10.1002/gamm.201310002/full, as well as the [Matlab code in Stefan Güttel's webpage](http://guettel.com/) which you can use to compute it.

All these algorithms are approximate (but in the end, what is *not* approximate on a computer?), so you may encounter trouble if the element that you wish to compute is many orders of magnitudes smaller than $\|\exp(A)\|$. In that case, you may be interested in a bound on the decay rate of off-diagonal entries in https://arxiv.org/abs/1501.07376.