We analyze the complexity of the Walk on Spheres algorithm for simulating Brownian Motion in a domain Omega in R^d. The algorithm, which was first proposed in the 1950s, produces samples from the hitting probability distribution of the Brownian Motion process on boundary of Omega within an error of epsilon. The algorithm is used as a building block for solving a variety of differential equations, including the Dirichlet Problem. The WoS algorithm simulates a BM starting at a point X_0=x in a given bounded domain Omega until it gets epsilon-close to the boundary of Omega. At every step, the algorithm measures the distance d_k from its current position X_k to the boundary of Omega and jumps a distance of d_k/2 in a uniformly random direction from X_k to obtain X_{k+1}. The algorithm terminates when it reaches X_n that is epsilon-close to the boundary of Omega. It is not hard to see that the algorithm requires at least Omega(\log 1/epsilon) steps to converge. Only partial results with respect to upper bounds existed. In 1959 M. Motoo established an O(\log 1/epsilon) bound on the running time for convex domains. The results were later generalized for a wider, but still very restricted, class of planar and 3-dimensional domains by G.A. Mikhailov (1979). In our earlier work (2007), we established an upper bound of O(\log^2 1/epsilon) on the rate of convergence of WoS for arbitrary planar domains. We introduce subharmonic energy functions to obtain very general upper bounds on the convergence of the algorithm. Special instances of the upper bounds yield the following results for bounded domains Omega: * if Omega is a planar domain with connected exterior, the WoS converges in O(\log 1/epsilon) steps; * if Omega is a domain in R^3 with connected exterior, the WoS converges in O(\log^2 1/epsilon) steps; * for d>2, if Omega is a domain in R^d, the WoS converges in O((1/epsilon)^{2-4/d}) steps; * for d>3 if Omega is a domain in R^d with connected exterior, the WoS converges in O((1/epsilon)^{2-4/(d-1)}) steps; * for any d if Omega is a domain in R^d bounded by a smooth surface, the WoS converges in O(\log 1/epsilon) steps. We also demonstrate that the bounds are tight, i.e. we construct a domain from each class for which the upper bound is exact. Our results give the optimal upper bound of O(\log1/epsilon) in many cases for which only a bound polynomial in 1/epsilon was previously known. 30 minutes would be fine for the talk. I prefer not to talk on Monday due to jetlag.