We develop a new methodology for the efficient computation of epidemic final size distributions for a broad class of Markovian models. We exploit a particular representation of the stochastic epidemic process to derive a method which is both computationally efficient and numerically stable. The algorithms we present are also physically transparent and so allow us to extend this method from the basic SIR model to a model with a phase-type infectious period and another with waning immunity. The underlying theory is applicable to many Markovian models where we wish to efficiently calculate hitting probabilities.