We present a theoretical investigation of the elastic interactions in a heteroepitaxial bilayer consisting of a (001) tetragonal PbZr (x) Ti1-x O-3 and (001) rhombohedral PbZr1-x Ti (x) O-3 on a thick (001) passive substrate. Analytical expressions for the elastic interaction energies between the layers and the resultant ferroelastic twin formation have been derived as a function of the lattice misfit strain (between layers and the substrate), composition of the ferroelectric and thickness. It is found that the elastic coupling between the tetragonal and rhombohedral layers leads to the equilibrium domain fraction in the tetragonal layer several time larger than that in single-layer films of similar thickness. Most critically, the model finds a significant change in the ferroelastic domain volume fraction in the presence of an applied electric field and hence enhanced piezoelectric properties compared to single-layered epitaxial PZT thin films.