Based on characteristic functions of variants, we developed an unconventional phase field modeling for investigating domains formation and evolution in tetragonal ferroelectrics. In order to develop this computational approach, we constructed the anisotropy energy of tetragonal variants, which is used instead of Landau-Devonshire potential in the conventional phase field method, resulting in that much fewer parameters are needed for simulations. This approach is advantageous in simulations of emerging ferroelectric materials. We employ it to study the formation and evolution of domains in tetragonal barium titanate single crystal, as well as the nonlinear behaviors under cyclical stress and electric field loading. A multi-rank laminated ferroelectric domain pattern, 90° domain switching accompanied by polarization rotation, and 180° domain switching accompanied by move of domain wall are predicted. It is found that the speed of 90° domain switching is slower than that of 180° domain switching, due to both polarization and transformation strain changed in 90° domain switching. It also suggests that large strain actuation can be generated in single crystal ferroelectrics via combined electromechanical loading inducing 90° domain switching. The good agreement between simulation results and experimental measurements is observed. Based on characteristic functions of variants, we developed an unconventional phase field modeling for investigating domains formation and evolution in tetragonal ferroelectrics. We order the anisotropy energy of tetragonal variants, which is used instead of Landau-Devonshire potential in the conventional phase field method, resulting in that that much fewer parameters are needed for simulations. This approach is advantageous in simulations of emerging ferroelectric materials. We employ it to study the formation and evolution of domains in tetragonal barium titanate single crystal, as well as the nonlinear behaviors under cyclical stress and electric field loading. A multi-rank laminated ferroelectric domain pattern, 90 ° domain switching accompanied by polarization rotation, and 180 ° domain switching accompanied by move of domain wall are predicted. It is found that the speed of 90 ° domain switching is slower than that of 180 ° domai n switching, due to both polarization and transformation strain changed in 90 ° domain switching. It also suggests that large strain actuation can be generated in single crystal ferroelectrics via combined crystal ferroelectrics via combined electromechanical loading inducing 90 ° domain switching. is observed.