Despite the undeniable success of the standard model of particle physics (SM), there are some phenomena that the SM cannot explain. These phenomena indicate that the SM has to be modified. One of the possible ways to extend the SM is to introduce heavy neutral leptons (HNLs). To search for HNLs in intensity Frontier experiments, one has to consider HNL production both in two-body and three-body decays of some mesons. We verified the possibility of using the parton level PYTHIA default matrix elements (without the form-factor formalism) to calculate HNL production in three-body semileptonic decays of B and D mesons in the experimentally interesting mass range of the produced HNLs. We conclude that this approach is quite suitable for the estimation of the sensitivity region for HNLs in the intensity Frontier experiments, provided one uses suitable parton level PYTHIA default matrix elements. Our study was driven by the usage of such an approximation by the SHiP collaboration. We conclude that in this case the parton level PYTHIA default matrix elements could have been chosen more appropriately.