We study a numerical method for solving a system of Volterra-renewal integral equations with space fluxes, that represents the Chapman-Kolmogorov equation for a class of piecewise deterministic stochastic processes. The solution of this equation is related to the time dependent distribution function of the stochastic process and it is a non-negative and non-decreasing function of the space. Based on the Bernstein polynomials, we build up and prove a non-negative and non-decreasing numerical method to solve that equation, with quadratic convergence order in space.