For non-linear solution of the Boltzmann equation, the cumulant moment method has been studied. To apply the method to the normal shock wave problem, we restricted ourselves to the monatomic Maxwell molecular gases. The method is based on the iterative approach developed by Maxwell-Ikenberry-Truesdell (MIT). The original MIT approach employs the equilibrium distribution function for the initial values in beginning the iteration. In the present work, we use the Mott-Smith bimodal distribution function to calculate the initial values and follow the MIT iteration procedure. Calculations have been carried out up to the second iteration for the profiles of density, temperature, stress, heat flux, and shock thickness of strong shocks, including the weak shock thickness of Mach range less than 1.4. The first iteration gives a simple analytic expression for the shock profile, and the weak shock thickness limiting law which is in exact accord with the Navier-Stokes theory. The second iteration shows that the calculated strong shock profiles are consistent with the Monte Carlo values quantitatively.