The Bessel beam belongs to a typical class of non-diffractive optical fields that are characterized by their invariant transverse profiles with the beam propagation. The extended and uniformed intensity profile in the axial direction is of great interest in many applications. However, ideal Bessel beams only rigorously exist in theory; the Bessel beams generated in the experiment are always quasi-Bessel beams with finite focal extensions and varying intensity profiles along the propagation axis. The ability to shape the on-axis intensity profile to meet specific needs is essential for many applications. Here, we demonstrate an iterative optimization based approach to engineer the on-axis intensity of Bessel beams through design and fine-tune processes. Starting with a standard axicon phase mask, the design process uses the computed on-axis beam profile as a feedback in the iterative optimization process, which searches for the optimal radial phase distribution that can generate a so-called generalized Bessel beam with the desired on-axis intensity profile. The fine-tune process repeats the optimization processing by using the adjusted target on-axis profile according to the measured one. Our proposed method has been demonstrated in engineering several quasi-Bessel beams with customized on-axis profiles. The high accuracy and high energy throughput merit its use in many applications. This method is also suitable to engineer higher-order Bessel beams by adding appropriate vortex phases into the designed phase mask. © 2018 IOP Publishing Ltd.