  We consider Monte Carlo methods for simulating solutions to the analogue of
the Dirichlet boundary-value problem in which the Laplacian is replaced by the
fractional Laplacian and boundary conditions are replaced by conditions on the
exterior of the domain. Specifically, we consider the analogue of the so-called
`walk-on-spheres` algorithm. In the diffusive setting, this entails sampling
the path of Brownian motion as it uniformly exits a sequence of spheres
maximally inscribed in the domain. As this algorithm would otherwise never end,
it is truncated when the `walk-on-spheres` comes within epsilon > 0 of the
boundary. In the setting of the fractional Laplacian, the role of Brownian
motion is replaced by an isotropic alpha-stable process with alpha in (0, 2). A
significant difference to the Brownian setting is that the stable processes
will exit spheres by a jump rather than hitting their boundary. This difference
ensures that disconnected domains may be considered and that, unlike the
diffusive setting, the algorithm ends after an almost surely finite number of
steps.
