The semi-implicit Euler-Maruyama (EM) method is investigated to approximate a class of time-changed stochastic differential equations, whose drift coefficient can grow super-linearly and diffusion coefficient obeys the global Lipschitz condition. The strong convergence of the semi-implicit EM is proved and the convergence rate is discussed. When the Bernstein function of the inverse subordinator (time-change) is regularly varying at zero, we establish the mean square polynomial stability of the underlying equations. In addition, the numerical method is proved to be able to preserve such an asymptotic property. Numerical simulations are presented to demonstrate the theoretical results.