This paper models the domain dynamics in a ferroelastic epilayer within the time-dependent Ginzburg-Landau (TDGL) framework. Constrained on a paraelastic substrate of square symmetry, the epilayer has rectangular symmetry, and forms domains of two variants. The domain wall energy drives the domains to coarsen. The spontaneous strains induce an elastic field, which drives the domains to refine. The competition between coarsening and refining selects an equilibrium domain size. We model the epilayer-substrate as a nonequilibrium thermodynamic system, evolving by the changes in the elastic displacements and the order parameters. The free energy consists of two parts: the bulk elastic energy, and the excess surface energy. The elastic energy density is taken to be quadratic in the strains. The surface energy density is expanded into a polynomial of the order parameters, the gradients of the order parameters, and the strains. In this expansion, the surface stress is taken to be quadratic in the order parameters. The evolution equations are derived from the free energy variation with respect to the order parameters. The elastic field is determined by superposing the Cerruti solution. Examples of computer simulation are presented.