The one dimensional transient heat conduction equation was numerically modeled through matrix diagonalization and three time-discretization schemes. The discrete methods were first-order backward, second-order backward, and implicit finite difference schemes. All simulations used the central difference formula in the space dimension. Two relevant physical systems were considered: a uniformly conducting slab and a melting block of ice. The latter lead to a moving boundary system, or Stefan problem. The multiphysics of melting was numerically modeled through alternating updates of temperature and melt front profiles. Iterative simulations were run with regularly refined discretization meshes in both systems. In the case of the conducting slab, temperature at a fixed point in space and time was considered. For the Stefan problem, the melt front movement after a set time was the physical solution of interest. The accuracy of the convergent results was increased using Richardson acceleration and the Wynn's epsilon algorithm. Accuracy was improved for the moving boundary problem as well, but to a significantly lesser degree. The relative errors improved by five and two orders of magnitude for the conducting block and melting ice simulations, respectively. These relative errors were used to determine that matrix diagonalization is the most accurate numerical solution among the four considered. In both simulation convergence and acceleration potential, matrix diagonalization was superior to the implicit and explicit discretization solutions. However, matrix diagonalization required significantly more computational time. With the enhancement of convergence acceleration, the finite difference schemes obtained similar relative errors to the diagonalization model. This demonstrated the value of convergence acceleration in the classic dilemma for every programmer. There is always a balance struck between model sophistication, accuracy, and computational time. Convergence acceleration allows for a simpler numerical model to achieve comparable accuracy, and in less time than that required for sophisticated numerical models. The numerical models were also compared for stability through parameters that arose in each simulation. These parameters were the Courant-Friederichs-Lewy (CFL) condition and diagonalized eigenvalues. Though diagonalization was found to be the most accurate, it was determined that the backwards finite difference solutions are the easiest to evaluate for stability. In these solution methods, the CFL value allows the stability to be determined prior to running the simulation.