A system is differentially flat if it is Lie-Bäcklund (L-B) equivalent to a free dynamical system that has dimensions equal to that of the input of the original system. Utilizing this equivalence, the problem of nonlinear model predictive control of a flat system can be reduced to a lower dimensional nonlinear programming problem with respect to the flat outputs. In this work, a novel computational method based on Haar wavelets in the time-domain for solving the resulting nonlinear programming problem is developed to obtain an approximation of the optimal flat output trajectory. The Haar wavelet integral operational matrix is utilized to transform the nonlinear programming problem to a finite dimensional nonlinear optimization problem. The proposed approach makes use of flatness as a structural property of nonlinear systems and the convenient mathematical properties of Haar wavelets to develop an efficient computational algorithm for nonlinear model predictive control of differentially flat systems. Further improvement on computational efficiency is achieved by providing solutions with multiple resolutions (e.g., obtaining high resolution solutions only for the near future, but allowing coarse approximation for the later stage in the prediction horizon).