Forward-time, centered space evalaution of the heat equation: numerical stability and unique solution

Question

I have a script of code which models a planetesimal that is accreted into a planetary atmosphere. In the code, I include the physics of frictional ablation and thermal ablation. Frictional ablation is governed by the coupled evolution of 5 equations and thermal ablation is governed by sublimatting off layers of the planetesimal once the ambient temperature in the atmosphere is greater than the planetesimal material's melting point. I am having a longstanding problem, however, with which I seek sincere assistance. I know the frictional ablation calculation is fine, but once thermal ablation comes in there are problems. I am calculating thermal ablation by using the forward-time, centered-space finite-difference method. Hence, this requires values for the time and spatial steps, dt and dr, respectively. However, for different values of these finite-differences, I get significantly different solutions for my thermal ablation profile in the output (figure 114 in the code). This is obviously a problem since I require an approximative convergence to a unique solution. I wonder if there are expert numerical analysts that may be willing to offer their knowledge into a solution. I want to include my python code and a text file for reference. To proceed, does anyone know how I can create a link online to share these with anyone?