A computer program has been developed which enables us to calculate the non-Fickian diffusion equations based on the Caputo fractional derivative. In recent years, the geological disposal method is, among to others, one of the most promising ways to dispose high-level radioactive wastes. It is therefore necessary to develop models for predicting solute transport in subsurface rock masses in order to assess and ensure the reliability for the geological disposal method. The conventional mathematical model of solute transport in rock masses is based on the Fick's law. However, because of existing complex fractures in the rock matrix, the conventional model tends to predict smaller solute travel distance than that in the actual transport process. The deficiency in the conventional model calls for the necessity of a more quantitative model which is applicable to the fractured porous media. In contrast, the non-Fickian diffusion model can provide realistic representation of actual fluid flow in the heterogeneous media. In the non-Fickian diffusion, the diffusion equation is described by fractional-in-space derivative of order α, which may vary from 0 to 1. In this study, we take advantage of the Caputo fractional derivative equation in order to develop a numerical method for analyzing the non-Fickian diffusion. We provide a numerical solution of the fractional derivative equation using implicit-finite difference method. The numerical result obtained for one dimensional diffusion equation using the computer program was shown to be in a good agreement with the analytical solution. The numerical method developed in this study was also verified with the result of laboratory experiment, in which a thermoluminescence technique was employed to trace the solute transport front in a granite sample. It has been shown that the numerical method based on the non-Fickian diffusion equation provides a better characterization of the experimental result, compared with the conventional diffusion equation.