The nonlinearity plays an important role in modulation of wave propagation characteristics, such as the adjustment of band structures by media nonlinearities. There was a modified incremental harmonic balance (IHB) method developed to hand single- and two-degree-of-freedom (2-DoF) wave propagation problems with nonlinearities. When the nonlinearity is strong, high-order basis functions are used to achieve accurate solutions, which requires complex formula derivation and high computational effort, especially for a general multi-DoF wave problem. To handle the issue, an improved IHB method is developed for solving multi-DoF wave problems with strong nonlinearities. In the method, construction of analytical Jacobian matrices of delay differential equations that are transformed from original wave problems is formulated for any DoFs and expansion orders, by which no manual derivation and numerical integration are needed. Thus, computational effort can be significantly reduced, and convergence analysis of solutions obtained from the IHB method with different expansion orders can be conducted to determine the minimum order of stable solutions, which has not been done for wave problems with strong nonlinearities. Moreover, the algorithm to solve nonlinear wave problems with given frequencies and unknown wave vectors is newly proposed in this work based on the improved IHB method, which is can be a more efficient approach if external excitations are given. In this work, one-dimensional diatomic lattices with Hertz contact law and cubic spring models are studied as examples, where dispersion curves and bandgap properties are generated. Results show that it is necessary to use high-order Fourier functions to obtain accurate steady-state wave responses when strong nonlinearities exist, which was not explicitly recognized in past works without convergence analysis for wave problems. Results also show that the convergence rate of the method with use of a fixed frequency is faster than that of a fixed wave vector.