We prove and implement stochastic solution (or Feynman-Kac) formulas for boundary value problems involving the spectral fractional Laplacian with nonzero Dirichlet boundary condition. The main tools used in the proofs are the abstract Cauchy problem for Feller semigroups together with Balakrishnan's theory of fractional powers. We show the the spectral fractional Laplacian with nonzero Dirichlet boundary conditions is the generator of an appropriate Feller semigroup for subordinate stopped Brownian motion, and obtain a stochastic solution formula for the fractional heat equation in a bounded domain. We then obtain precise regularity and steady-state convergence properties of the parabolic problem using the eigenfunction expansion of the classical solution, which leads to estimates for the survival probability of subordinate stopped Brownian motion. These results allow us to take the long-time limit in the parabolic formula to establish a stochastic solution formula for the Dirichlet boundary value problem. These stochastic solution formulas for the spectral fractional Laplacian with nonzero boundary conditions are novel, and allow for efficient, embarrassingly parallel local solution of the above boundary value problems. We discuss the discretization of these formulas, and verify them in dimensions two and three with benchmark examples. We study the effect of the number of path samples and the time step size for path discretization on the accuracy of the solution.