We consider the path approximation of Bessel processes and develop a new and efficient algorithm. This study is based on a recent work by the authors, on the path approximation of the Brownian motion, and on the construction of specific own techniques. It is part of the family of the so-called $\varepsilon$-strong approximations. More precisely, our approach constructs jointly the sequences of exit times and corresponding exit positions of some well-chosen domains, the construction of these domains being an important step. Based on this procedure, we emphasize an algorithm which is easy to implement. Moreover, we can develop the method for any dimension. We treat separately the integer dimension case and the non integer framework, each situation requiring appropriate techniques. In particular, for both situations, we show the convergence of the scheme and provide the control of the efficiency with respect to the small parameter $\varepsilon$. We expand the theoretical part by a series of numerical developments.