The purpose of this paper is to design a fully discrete hybridized discon-tinuous Galerkin (HDG) method for solving a system of two-dimensional (2D) coupled Burgers equations over a specified spatial domain. The semi-discrete HDG method is designed for a nonlinear variational formulation on the spatial domain. By exploiting broken Sobolev approximation spaces in the HDG scheme, numerical fluxes are defined properly. It is shown that the proposed method is stable under specific mild conditions on the stabi-lization parameters to solve a well-posed (in the sense of energy method) 2D coupled Burgers equations, which is imposed by Dirichlet boundary conditions. The fully discrete HDG scheme is designed by exploiting the Crank–Nicolson method for time discretization. Also, the Newton–Raphson method that has the order of at least two is nominated for solving the obtained nonlinear system of coupled Burgers equations over the rect-angular domain. To reduce the complexity of the proposed method and the size of the linear system, we exploit the Schur complement idea. Numerical results declare that the best possible rates of convergence are achieved for approximate solutions of the 2D coupled Burgers equations and their first-order derivatives. Moreover, the proposed HDG method is examined for two other types of systems, that is, a system with high Reynolds numbers and a system with an unavailable exact solution. The acceptable results of examples show the flexibility of the proposed method in solving various problems
