We provide a novel expression of the scale function for a Lévy process with negative phase-type jumps. It is in terms of a certain transition rate matrix which is explicit up to a single positive number. A monotone iterative scheme for the calculation of the latter is presented and it is shown that the error decays exponentially fast. Our numerical examples suggest that this algorithm allows us to employ phase-type distributions with hundreds of phases, which is problematic when using the known formula for the scale function in terms of roots. Extensions to other distributions, such as matrix-exponential and infinite-dimensional phase-type, can be anticipated.
