A biochemical network can be simulated by a set of ordinary differential equations (ODE) under well-stirred reactor conditions, for large numbers of molecules, and frequent reactions. This is no longer a robust representation when some molecular species are in small numbers and reactions changing them are infrequent. In this case, discrete stochastic events trigger changes of the smooth deterministic dynamics of the biochemical network. Piecewise-deterministic Markov processes (PDMP) are well adapted for describing such situations. Although PDMP models are now well established in biology, these models remain computationally challenging. Previously we have introduced the push-forward method to compute how the probability measure is spread by the deterministic ODE flow of PDMPs, through the use of analytic expressions of the corresponding semigroup. In this paper we provide a more general simulation algorithm that works also for non-integrable systems. The method can be used for biochemical simulations with applications in fundamental biology, biotechnology and biocomputing. This work is an extended version of the work presented at the conference CMSB2019.