One very practical method of finding the $q$'s you are looking for is presented here:
Kovács, Attila, and Norbert Tihanyi. "Efficient computing of n-dimensional simultaneous Diophantine approximation problems." Acta Universitatis Sapientiae, Informatica 5.1 (2013): 16-34.
In what has to be one of the coolest applications of simultaneous Diophantine approximation ever, the authors use their algorithm to find many (but not all) values of the Riemann-Siegel Z function where Z(t) >500. See: http://riemann-siegel.com/index.php for more details.
Another choice is the method described in:
Bosma, Wieb, and Ionica Smeets. "Finding simultaneous Diophantine approximations with prescribed quality." The Open Book Series 1.1 (2013): 167-185.
From the abstract: "We give an algorithm that finds a sequence of approximations with Dirichlet coefficients bounded by a constant only depending on the dimension. The algorithm uses LLL lattice basis reduction."
You may find this paper helpful:
Lagarias, J. C. "Geodesic multidimensional continued fractions." Proceedings of the London Mathematical Society 3.3 (1994): 464-488.
Doug Hensley references the Lagarias algorithm in this paper (book chapter?):
If you are not too concerned with the quality of the approximations this book presents many algorithms that work (w/o LLL) and may be useful:
Brentjes, Arne Johan. "Multi-dimensional continued fraction algorithms." MC Tracts 145 (1981): 1-183.
This book may still be available for purchase through ntis.gov and is well worth the cost.
A complete answer to your question depends on how the $q$'s you are looking for compare to the given $\epsilon$ and the desired computational complexity of the process for finding them. In the extremes, I don't know if there is a completely satisfying answer.