Abstract: The Method of Characteristics (MOC) is capable to accurately solve the neutron transport equation with arbitrary geometry. However, the MOC suffers from some drawbacks: slow convergence and time consuming. Based on the spatial domain decomposition and the ray parallelization, the parallel 2D MOC algorithm was implemented with MPI+OepnMP/CUDA programming model to leverage the computing power of Central Processing Unit-Graphics Processing Unit (CPU-GPU) heterogeneous high-performance computing systems. In addition, a dynamic workload partitioning scheme was proposed to efficiently take advantage of all the CPU and GPU resources. The workload is appropriately assigned to the CPU and GPU according to their computational capabilities, and all CPUs and GPUs perform the calculation concurrently. The numerical results demonstrate that the parallel algorithm maintains the desired accuracy. Meanwhile, the dynamic workload portioning scheme can provide the optimal workload partition based on the runtime performance. As a result, about 14% improvement is observed in the overall performance compared with the MPI+CUDA parallelization when the CPU-GPU heterogeneous computation is performed on 5 heterogeneous nodes (including 20 GPUs).