First exit times of stochastic processes are fundamental in many applications. In mathematical finance, they are used to quantify default risk in path-dependent derivatives; in neuroscience, they describe interspike interval distributions. Diffusion processes, as solutions of stochastic differential equations, form a central class of models, making the accurate approximation of their exit times a problem of broad interest.

We consider the multidimensional setting and study the numerical approximation of the first exit time \(\tau_{\mathcal D}\) of a \(d\)-dimensional diffusion process \((X_t)_{t \ge 0}\) from a bounded, regular domain \(\mathcal D\). The process satisfies \[ dX_t = \nabla \mathcal{U}(X_t,t),dt + dB_t, \qquad X_0 \in \mathcal D, \] where \((B_t)\) is a \(d\)-dimensional Brownian motion and the drift term may depend on both space and time. Our objective is to design an efficient alternative to the classical Euler scheme, which requires small time steps to ensure accuracy near the boundary.

In the Brownian case, the Random Walk on Spheres (WOS) algorithm exploits isotropy to perform large spatial jumps, leading to a mean number of steps proportional to \(|\log(\varepsilon)|\), where \(\varepsilon\) is the boundary layer parameter. Extensions based on spheroids allow the joint approximation of exit position and exit time.

We generalize this approach to multidimensional diffusion processes with drift. The proposed method relies on an acceptance–rejection procedure applied to random walk trajectories and introduces truncated spheroids to account for nonzero drift. This construction preserves the efficiency of large spatial displacements while incorporating the effect of the drift term. The performance of the algorithm is supported by theoretical results and illustrated through numerical experiments.