Bilayer plates are slender structures made of two thin layers of different materials. They react to environmental stimuli and undergo large bending deformations with relatively small actuation. The reduced model is a constrained minimization problem for the second fundamental form, with a given spontaneous curvature that encodes material properties, subject to an isometry constraint. We design a local discontinuous Galerkin (LDG) method, which imposes a relaxed discrete isometry constraint and controls deformation gradients at barycenters of elements. We prove Gamma-convergence of LDG, design a fully practical gradient flow, which gives rise to a linear scheme at every step, and show energy stability and control of the isometry defect. We extend the Gamma-convergence analysis to piecewise quadratic creases. We also illustrate the performance of the LDG method with several insightful simulations of large deformations, one including a curved crease.