  We introduce a general algorithm for the computation of the scale functions
of a spectrally negative L\'evy process $X$, based on a natural weak
approximation of $X$ via upwards skip-free continuous-time Markov chains with
stationary independent increments. The algorithm consists of evaluating a
finite linear recursion with its (nonnegative) coefficients given explicitly in
terms of the L\'evy triplet of $X$. Thus it is easy to implement and
numerically stable. Our main result establishes sharp rates of convergence of
this algorithm providing an explicit link between the semimartingale
characteristics of $X$ and its scale functions, not unlike the one-dimensional
It\^o diffusion setting, where scale functions are expressed in terms of
certain integrals of the coefficients of the governing SDE.
