This paper explores and implements high-order numerical schemes for integrating linear parabolic partial differential equations with piece-wise smooth boundary data. The high-order Monte-Carlo methods we present give extremely accurate approximations in computation times that we believe are comparable with much less accurate finite difference and basic Monte-Carlo schemes. A key step in these algorithms seems to be that the order of the approximation is tuned to the accuracy one requires. A considerable improvement in efficiency can be attained by using ultra high-order cubature formulae. Lyons and Victoir (“Cubature on Wiener Space, Proc. R. Soc. Lond. A 460, 169–198”) give a degree 5 approximation of Brownian motion. We extend this cubature to degrees 9 and 11 in 1-dimensional space-time. The benefits are immediately apparent.