A new scheme based on the homotopy analysis method (HAM) is developed for calculating the nonlinear normal modes (NNMs) of multi degrees-of-freedom (MDOF) oscillatory systems with quadratic and cubic nonlinearities. The NNMs in the presence of internal resonances can also be computed by the proposed method. The method starts by approximating the solution at the zeroth-order, using some few harmonics, and proceeds to higher orders to improve the approximation by automatically including higher harmonics. The capabilities and limitations of the method are thoroughly investigated by applying them to three nonlinear systems with different nonlinear behaviors. These include a two degrees-of-freedom (2DOF) system with cubic nonlinearities and one-to-three internal resonance that occurs on nonlinear frequencies at high amplitudes, a 2DOF system with quadratic and cubic nonlinearities having one-to-two internal resonance, and the discretized equations of motion of a cylindrical shell. The later one has internal resonance of one-to-one. Moreover, it has the symmetry property and its DOFs may oscillate with phase difference of 90 deg, leading to the traveling wave mode. In most cases, the estimated backbone curves are compared by the numerical solutions obtained by continuation of periodic orbits. The method is found to be accurate for reasonably high amplitude vibration especially when only cubic nonlinearities are present.