In many modern applications of queueing theory, the classical assumption of exponentially decaying service distributions does not apply. In particular, Internet and insurance risk problems may involve heavy-tailed distributions. A difficulty with heavy-tailed distributions is that they may not have closed-form, analytic Laplace transforms. This makes numerical methods, which use the Laplace transform, challenging. In this paper, we develop a method for approximating Laplace transforms. Using the approximation, we give algorithms to compute the steady state probability distribution of the waiting time of an M/G/1 queue to a desired accuracy. We give several numerical examples, and we validate the approximation with known results where possible or with simulations otherwise. We also give convergence proofs for the methods.