In electron probe microanalysis or scanning electron microscopy, the Monte Carlo method is widely used for modeling electron transport within specimens and calculating X-ray spectra. For an accurate simulation, the calculation of secondary fluorescence (SF) is necessary, especially for samples with complex geometries. In this study, we developed a program, using a hybrid model that combines the Monte Carlo simulation with an analytical model, to perform SF correction for three-dimensional (3D) heterogeneous materials. The Monte Carlo simulation is performed using MC X-ray, a Monte Carlo program, to obtain the 3D primary X-ray distribution, which becomes the input of the analytical model. The voxel-based calculation of MC X-ray enables the model to be applicable to arbitrary samples. We demonstrate the derivation of the analytical model in detail and present the 3D X-ray distributions for both primary and secondary fluorescence to illustrate the capability of our program. Examples for non-diffusion couples and spherical inclusions inside matrices are shown. The results of our program are compared with experimental data from references and with results from other Monte Carlo codes. They are found to be in good agreement.