We propose a second order discretization for backward stochastic differential equations (BSDEs) with possibly nonsmooth boundary data. When implemented, the discretization method requires essentially the same computational effort with the Euler scheme for BSDEs of Bouchard and Touzi [Stochastic Process. Appl. 111 (2004) 175–206] and Zhang [Ann. Appl. Probab. 14 (2004) 459–488]. However, it enjoys a second order asymptotic rate of convergence, provided that the coefficients of the equation are sufficiently smooth. In the second part of the paper, we combine this discretization with higher order cubature formulas on Wiener space to produce a fully implementable second order scheme.