A new fully numerical method is presented which employs multiple Poincaré sections to find quasi-periodic orbits. The main advantages of this method are the small overhead cost of programming and very fast execution times, robust behavior near chaotic regions that leads to full convergence for given family of quasi-periodic orbits and the minimal memory required to store these orbits. This method reduces the calculation of the search for the two-dimensional invariant torus to a search for the closed orbits, which are the intersection of the invariant torus with the Poincaré sections. Truncated Fourier series are employed to represent these closed orbits. The flow of the differential equation on the invariant torus is reduced to maps between the consecutive Poincaré maps. A Newton iteration scheme makes use of the invariancy of the circles of the maps on these Poincaré sections in order to find the Fourier coefficient that define the circles to any given accuracy. A continuation procedure that uses the incremental behavior of the Fourier coefficients between close quasi-periodic orbits is utilized to extend the results from a single orbit to a family of orbits. Quasi-Halo and Lissajous families of the Sun-Earth Restricted Three-Body Problem (RTBP) around the L1 and L2 libration points are obtained via this method. Results are compared with the existing literature.