The analytical solutions of the one- and two-phase Stefan problems are found in the form of series containing linear combinations of the integral error functions which satisfy a priori the heat equation. The unknown coefficients are defined from the initial and boundary conditions by the comparison of the like power terms of the series using the Faa di Bruno formula. The convergence of the series for the temperature and for the free boundary is proved. The approximate solution is found using the replacement of series by the corresponding finite sums and the collocation method. The presented test examples confirm a good approximation of such approach. This method is applied for the solution of the Stefan problem describing the dynamics of the interaction of the electrical arc with electrodes and corresponding erosion.