The bending of bilayer plates is a mechanism that allows for large deformations via small externally induced lattice mismatches of the underlying materials. Its mathematical modeling, discussed herein, consists of a nonlinear fourthorder problem with a pointwise isometry constraint. A discretization based on Kirchhoff quadrilaterals is devised and its -convergence is proved. An iterative method that decreases the energy is proposed, and its convergence to stationary configurations is investigated. Its performance, as well as reduced model capabilities, are explored via several insightful numerical experiments involving large (geometrically nonlinear) deformations.