We develop a discretization method for solving the minimal energy configuration of bilayer plates based on a mathematical model developed in Schmidt (2007), Bartels et al. (2017). Our discretization method employs C1-spline functions. A highlight of the method involves a trick to handle the nonlinear isometry constraint in such a way that not only numerical integration becomes unnecessary, but also that the final optimization problems are in the form of degree 4 polynomial optimization problems (POP). We develop two different versions of the method, one resulted in a constrained degree 4 POP involving a small tolerance ɛ, another resulted in an unconstrained degree 4 POP involving a large penalty parameter μ. We develop a mathematical analysis, based on the direct method and techniques in Γ-convergence, to show how ɛ and μ can be chosen according to the grid size so that the minimizers of the discrete problems converge to that of the continuum variational problem as the grid size goes to zero. We corroborate the theory through a series of computational experiments, and also report an unexpected finding related to the asymmetry of the discretized problems.