A phase-space distribution function of the steady state in galaxy models that admits regular orbits overall in the phase-space can be represented by a function of three action variables. This type of distribution function in Galactic models is often constructed theoretically for comparison of the Galactic models with observational data as a test of the models. On the other hand, observations give Cartesian phase-space coordinates of stars. Therefore it is necessary to relate action variables and Cartesian coordinates in investigating whether the distribution function constructed in galaxy models can explain observational data. Generating functions are very useful in practice for this purpose, because calculations of relations between action variables and Cartesian coordinates by generating functions do not require a lot of computational time or computer memory in comparison with direct numerical integration calculations of stellar orbits. Here, we propose a new method called a torus-fitting method, by which a generating function is derived numerically for models of the Galactic potential in which almost all orbits are regular. We confirmed the torus-fitting method can be applied to major orbit families (box and loop orbits) in some two-dimensional potentials. Furthermore, the torus-fitting method is still applicable to resonant orbit families, besides major orbit families. Hence the torus-fitting method is useful for analyzing real Galactic systems in which a lot of resonant orbit families might exist.