In this paper we present an algorithm of computing two-dimensional (2D) stable and unstable manifolds of hyperbolic fixed points of nonlinear maps. The 2D manifold is computed by covering it with orbits of one-dimensional (1D) sub-manifolds. A generalized Foliation condition is proposed to measure the growth of 1D sub-manifolds and eventually control the growth of the 2D manifold along the orbits of 1D sub-manifolds in different directions. At the same time, a procedure for inserting 1D sub-manifolds between adjacent sub-manifolds is presented. The recursive procedure resolves the insertion of new mesh point, the searching for the image (or pre-image), and the computation of the 1D sub-manifolds following the new mesh point tactfully, which does not require the 1D sub-manifolds to be computed from the initial circle and avoids the over assembling of mesh points. The performance of the algorithm is demonstrated with hyper chaotic three-dimensional (3D) Hnon map and Lorenz system.