In this work, a new solution is developed for the problem of contaminant transport in a single fracture-matrix system, where the first-order reaction rate constants are different in both fracture and matrix. It takes a form of convolution with three functions as a basis to consider different transport mechanisms separately. The statistical nature of the three functions, as well as the interpretation of the solution as a marginal probability distribution in the case of no first-order reactions, allows us to develop a simple Time Domain Random Walk (TDRW) algorithm to calculate the breakthrough curves at a given point of observation downstream the fracture. Compared with the existing versions of the TDRW algorithm, it is superior not only in the physical reasoning and statistical interpretations but also in its numerical implementations. In addition, the developed algorithm can not only be used to estimate the distribution profile of the contaminant concentration along the fracture but also the concentration within the matrix, since the analytical solution to contaminant concentration in the matrix also takes a convolution form of three functions. Also, the distribution profile of contaminant concentration within the matrix can readily be determined by the use of our TDRW algorithm. To validate the developed algorithm, three benchmark cases are considered for either nuclide or colloid transport through a fractured rock. The results show that TDRW algorithm is superior to the Gaussian quadrature solution, but similar to inverse Laplace transform solution, in computational expense when nearly identical results are obtained. However, the Monte Carlo nature of the TDRW algorithm implies that the accuracy of the computational result is related to the number of particles applied in the simulation, which might make the obtained results fluctuated.