In this investigation, a model of tumor growth with a free boundary is studied. In this model, the mutations of tumor cells, which lead to heterogeneous tumors, are considered. An uncountable set of equations is used to show the density of tumor cells. Each member of this set presents the dynamic of a specific cell by the inclusion of an independent variable y. So, all equations are presented by one equation dependent on y. The transformation of a cell from one type to another happens with probability \(P(y_1,y_2)\) . To include the transformations, consumption of drug and nutrient by cells, and the effects of them on cells, integral terms are added to the equations. Time fractional order equations form this model to include memories as well. In order to see the behavior and radius of the tumor under the effects of drug and nutrient, the dynamics of sensitive and resistant cells, and so on, we need to solve the problem. Since obtaining the exact solutions is impossible, we must solve the problem numerically. But without stability and convergence analysis the solutions are not reliable. The spectral method is used to solve the equations. For this, the solutions are considered linear combinations of polynomials satisfying the boundary conditions with coefficients dependent on time and y. Integral terms are approximated using Legendre-Gauss nodes and time derivatives (of order \(\alpha \) ) are estimated with error \(O(h^{2-\alpha })\) , where h is the step size. Finally, coefficients are obtained using the collocation method. Then, we have proven the convergence and the stability of the presented numerical method. Numerical examples are presented to justify the theoretical statements. Some figures are also presented to illustrate the effects of drug, nutrients, and the coefficients of the model on the tumor cells (including resistant and sensitive cells) and radius of the tumor.