The theory of multistate template-directed reversible copolymerization is developed by extending the method based on iterated function systems to matrices, taking into account the possibility of multiple activation states instead of a single one for the growth process. In this extended theory, the mean growth velocity is obtained with an iterated matrix function system and the probabilities of copolymer sequences are given by matrix products defined along the template. The theory allows us to understand the effects of template heterogeneity, which include a fractal distribution of local growth velocities far enough from equilibrium, and a regime of sublinear growth in time close to equilibrium.