A numerical method for mean curvature motion in bounded domains with nonlinear Neumann boundary conditions is proposed and analyzed. It consists of a semi-Lagrangian scheme in the main part of the domain as proposed by Carlini, Falcone and Ferretti, combined with a finite difference scheme in small layers near the boundary to cope with the boundary condition. The consistency and monotonicity properties of the new scheme are studied for nonstructured triangular meshes in dimension two. Details on the implementation are given. Numerical tests are presented
