  We provide a novel expression of the scale function for a L\'evy processes
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 to employ phase-type distributions with a hundred 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.
