A semilinear reaction-diffusion equation with multiple solutions is considered in a smooth two-dimensional domain. Its diffusion parameter epsilon(2) is arbitrarily small, which induces boundary layers. Constructing discrete sub-and super-solutions, we prove existence and investigate the accuracy of multiple discrete solutions on layer-adapted meshes of Bakhvalov and Shishkin types. It is shown that one gets second-order convergence (with, in the case of the Shishkin mesh, a logarithmic factor) in the discrete maximum norm, uniformly in epsilon for epsilon <= Ch. Here h > 0 is the maximum side length of mesh elements, while the number of mesh nodes does not exceed Ch(-2). Numerical experiments are performed to support the theoretical results.