Among the numerous applications of the Ginzburg-Landau theory to the analysis of significant reaction-diffusion systems, modeling of the behavior of promising polar dielectric materials should be especially highlighted. The paper is devoted to the theoretical and numerical analysis of the Landau-Khalatnikov model describing the dynamics of 2D domain pattern formation in ferroelectrics. The unique solvability of the initial-boundary value problem for the system of 2D cubic-quintic Landau-Khalatnikov equations is proved. The proof is based on the derivation of new a priori estimates for the solution of the system of nonlinear parabolic equations. A series of computational experiments are conducted to examine both spontaneous and polar-induced domain pattern formation in biaxial ferroelectrics. Finite-element simulations allow us to visualize different types of ferroelectric domain structures depending on the varying boundary conditions.