If RS(n,k) code is defined as Sum[i=0..n-1] Xi Zi^j = 0, j=0..m-1, m=n-k, where {Xi} are codewords and {Zi} are locators, then for any polynomial P(Z) of degree less than m we have Sum Xi P(Zi) = 0.

Now, let's erase m values at locations E={Ej}, j=0..m-1. For each Ej, we can construct polynomial Pj(Z) of degree m-1 that evaluates to 1 at Ej and to 0 at all other locations in E. From Sum Xi Pj(Zi) = 0 we get Xj=Sum[i!=j] Xi Pj(Zi) (note that the sum is actually over all known Xi). {Pj(Zi)} is in fact the decoding matrix. Computing it is O(nm)=O(km^2) (not O(km) as I put above).

I have that algorithm implemented in my ErrorZ project on GitHub: https://github.com/OlegMazurov/ErrorZ/blob/master/src/main/j... More advanced stuff, like decoding up to m-1 errors, can also be found there.