Distribution-free tolerance regions are hyper-rectangles in terms of the number of variables that include at least a pre-specified proportion of normal subjects with a confidence bounded from below by a prescribed probability. This paper has two main goals. The first is to propose an innovative method for constructing multidimensional tolerance regions that work well when the only assumption that can be made about the underlying distribution is that it is continuous. Although our proposal is, in essence, an extension of the Wilks-Wald method to higher dimensions, this research is far less immediate and straightforward than it has appeared to many authors, who, moreover, have not really used it in practical work. In particular, the order statistics dividing the regions are not fixed beforehand, but are determined by an optimization procedure. The second goal is to suggest a way of overcoming a problem of practical importance concerning the order in which the variables are included, which has remained unsolved since the introduction of the argument almost eighty years ago. This device facilitates the search for a good solution while avoiding being drawn away into the black hole of combinatorial computations. Simulation and applications to laboratory medicine data illustrate the advantages of using the method presented in this paper.