The manuscript presents a nonlinear mathematical model capable to investigate the nonlinear bending response of Bi-directional functionally graded (BDFG) plate resting on elastic foundations including large deformations in the sense of von Kármán. A unified higher order shear plate theory is used to implement the shear influence with parabolic distribution through thickness direction. The gradation of materials is portrayed by power law function through axial (x-direction) and thickness (z-direction). The governing equations consist of a set of four nonlinear coupled partial differential equations. The assumption that the material properties change in the x- direction complicates the governing equations since they have variable-coefficients. Winkler–Pasternak model is exploited to present the elastic foundations with normal and shear deformations. In this work, a novel incremental–iterative method is proposed, and some tricky computational issues are performed to ensure stable and fast convergence rates even with large incremental load steps. The differential/integral quadrature method (DIQM) is developed to numerically discretize the variable-coefficient governing equations on a rectangular plate with different boundary conditions. Parametric studies are presented to illustrate the impact of nonlinear strains, gradation indices, geometrical properties, and boundary conditions on the transverse displacement, normal and shear stresses.