The paper describes the iteration method for finding the eigenfunctions and eigenvalues of the system of two nonlinear Schrödinger equations, which describes the process of second harmonic generation by femtosecond pulse in media with the quadratic and cubic nonlinear response. Coefficients, which characterize the nonlinearities, depend on one of the coordinate. The discussed method allows to find soliton solutions of new form corresponding to the first and second eigenvalues for the wide range of the nonlinear coefficients values. For determination of the eigenfunctions of the third and higher order it is necessary to select the initial approximation in a special way.