Many structures can behave nonlinearly, exhibiting behavior that is not captured by linear vibration theory such as localization and frequency-energy dependence. The nonlinear normal mode (NNM) concept, developed over the last few decades, can be quite helpful in characterizing a structure’s nonlinear response. In the definition of interest, an NNM is a periodic solution to the conservative nonlinear equations of motion. Several approaches have been suggested for computing NNMs and some have been quite successful even for systems with hundreds of degrees of freedom. However, existing methods are still too expensive to employ on realistic nonlinear finite element models, especially when the Jacobian of the equations of motion is not available analytically. This work presents a new approach for numerically calculating nonlinear normal modes by combining force appropriation, numerical integration and continuation techniques. This method does not require gradients, is found to compute the NNMs accurately up to moderate response amplitudes, and could be readily extended to experimentally characterize nonlinear structures. The method is demonstrated on a nonlinear mass-spring-damper system, computing its NNMs up to a 35% shift in frequency. The results are compared with those from a gradient based algorithm and the relative merits of each method are discussed.