A semilinear reaction-diffusion two-point boundary value problem, whose second-order derivative is multiplied by a small positive parameter epsilon(2), is considered. It can have multiple solutions. The numerical computation of solutions having interior transition layers is analysed. It is demonstrated that the accurate computation of such solutions is exceptionally difficult. To address this difficulty, we propose an artificial-diffusion stabilization. For both standard and stabilised finite difference methods on suitable Shishkin meshes, we prove existence and investigate the accuracy of computed solutions by constructing discrete sub- and super-solutions. Convergence results are deduced that depend on the relative sizes of epsilon and N, where N is the number of mesh intervals. Numerical experiments are given in support of these theoretical results. Practical issues in using Newton's method to compute a discrete solution are discussed.