We discuss two independent methods of solution of a master equation whose biased jump transition rates account for long jumps of Lévy-stable type and admit a Boltzmannian (thermal) equilibrium to arise in the large time asymptotics of a probability density function ρ(x, t). Our main goal is to demonstrate a compatibility of a direct solution method (an explicit, albeit numerically assisted, integration of the master equation) with an indirect pathwise procedure, recently proposed in [Physica A 392, 3485, (2013)] as a valid tool for a dynamical analysis of non-Langevin jump-type processes. The path-wise method heavily relies on an accumulation of large sample path data, that are generated by means of a properly tailored Gillespie’s algorithm. Their statistical analysis in turn allows to infer the dynamics of ρ(x, t). However, no consistency check has been completed so far to demonstrate that both methods are fully compatible and indeed provide a solution of the same dynamical problem. Presently we remove this gap, with a focus on potential deficiencies (various cutoffs, including those upon the jump size) of approximations involved in simulation routines and solutions protocols.