Weak-sense numerical methods for (underdamped) Langevin dynamics in bounded domains are constructed and analysed, with both their finite time convergence and convergence to ergodic limits being proved.  First-order methods are based on an Euler-type scheme interlaced with collisions with the boundary. To achieve second order,  composition schemes are derived based on decomposition of the generator into collisional drift,  impulse, and stochastic momentum evolution.  In a deterministic setting, this approach would typically lead to first-order approximation, even in symmetric compositions, but we find that the stochastic method can provide second-order weak approximation with a single force evaluation, both at finite times and in the ergodic limit.   We provide theoretical and numerical justification for this observation using model problems. A number of numerical experiments will be presented. The talk is based on a recent joint work with Ben Leimkuhler (Edinburgh) and Akash Sharma (Gothenburg).