Here is a way to get started with the numerical method presented in the paper.

Use your favorite coding language (Python, Matlab, Fortran...) to implement the code in Section 3.
This code is the code for the simplest finite difference method for the heat equation (38)
with formulas (39,40,41) to update values of s and of c at the last grid point
and with common sense rules for adding or removing grid points.
This should give you data for Table 5 right away.
Here is my version of the code written in Mathematica.

The simplest way to test the accuracy of the method is to compare the numerical solution
to known asymptotic values of the classical solution at small times. See Section 5.
To do this you only need to make sure that your code works correctly with at least 25 digit precision.
Here is my adaptation that produces raw data for Table 7.
Here is my code to obtain Richardson extrapolations from the raw data and gives Tables 7,8.
If you continue the calculation till the final time 0.1974 you can get Table 6.
Table 6 is based on calculations with various values of parameters, but with the smallest dx=0.00025.

Here is my calculation with dx=0.000125 which confirms old values in Table 6 and at many times suggests an extra correct digit.
This notebook also shows 'lstc' which lists all the pairs { x , c(x,0.1974) }, which together with the final value of x0 = s(0.1974)
can be used to continue the calculation till extinction (done here).
Thus various hypotheses about the behavior of the solution near extinction can be tested fast.

The numerical method was tested on how well it reproduces exact solutions in Examples 1,2,3.
My codes for Example 1 and Example 3 rely on Mathematica's ability to handle large integers easily.
If one wants to use Fortran, the evaluation of exact solutions for Examples 1,3 needs to be done differently.
However, no special considerations are needed for evaluation of the exact solution in Example 2.

Some experiments:

• Extinction. When the interface is sealed at x=0 (c_x(0,t)=0), one would expect that the solution c(.,t) goes to 0
  as t goes to infinity whenever c(x,0) > 0. However some odd behavior may sometime occur for some c(x,0)...
• Stability. Direct solutions of moving boundary problems with different initial values c(x,0), but with the same boundary values c(0,t), tend to approach the same solution. Here is an example where one solution is the solution in Example 1, determined by s(t) = 2- cos(t). The other solution is a direct numerical solution determined with same boundary values c(0,t), but totally different initial c(x,0) = 0.02 for 0 < x < 0.5. The differences between moving boundary positions and concentrations approach 0 fast. The same should work for Stefan free boundary problems.
• Resuscitation. When concentration c(x,0) is near 0, what kind of periodic conditions at x=0 are required
  for the concentration to revive. What is the good relation between amplitude and frequency?

If you have questions email me at milan@math.msu.edu