  Bernal-stacked bilayer graphene (BLG) provides an ideal basis for
gate-controlled, and free of etching, electronic devices. Theoretical modeling
of realistic devices is an essential part of research, however, simulations of
large-scale BLG devices continue to be extremely challenging. Micrometer-sized
systems are predominantly beyond the reach of the commonly used atomistic
tight-binding method, while other numerical approaches based on the two
dimensional Dirac equation are not straightforward to conduct due to the
fermion doubling problem. Here we present an approach based on the continuum
model, unharmed by the fermion doubling. The discretization of the BLG
continuum Hamiltonian leads to an effective four-band model, with both valleys
built-in. We demonstrate its performance with realistic, large-scale systems,
and obtain results consistent with experiments and with the tight-binding
model, over a broad range of magnetic field.
