This paper provides a computational method for dislocation-precipitate interaction problems. The computational method is a combination of the parametric dislocation dynamics and the voxel-based finite element method, and has a potential to enable the simulation of interaction between dislocations and multiple precipitates. To reduce the computational time, a multi-level voxel element model is employed. The convergence behavior of numerical accuracy and the computational time of the proposed method are examined by solving a dislocation-precipitate interaction problem. The results show that the proposed method has a good convergence behavior, and the computational time can be drastically reduced by the use of the multi-level voxel element model. Finally, the interaction between dislocations and multiple precipitates is solved to demonstrate a potential of the proposed method with various average diameters and constant volume fraction of precipitate. As the result, the proposed method successfully captured the dependence of the critical resolved shear stress on the average precipitate diameter.