  Recent studies of the fragmentation functions using the Nambu--Jona-Lasinio
(NJL) - Jet model have been successful in describing the quark fragmentation
functions to pions and kaons. The NJL-Jet model employs the integral equation
approach to solve for the fragmentation functions in quark-cascade description
of the hadron emission process, where one assumes that the initial quark has
infinite momentum and emits an infinite number of hadrons. Here we introduce a
Monte Carlo (MC) simulation method to solve for the fragmentation functions,,
that allows us to relax the above mentioned approximations. We demonstrate that
the results of MC simulations closely reproduce the solutions of the integral
equations in the limit where a large number of hadrons are emitted in the quark
cascade. The MC approach provides a strong foundation for the further
development of the NJL-Jet model that might include many more hadronic emission
channels with decays of the possible produced resonances, as well as inclusion
of the transverse momentum dependence (TMD), all of which are of considerable
importance to the experimental studies of the transverse structure of hadrons.
