Thermal diffusion of dopants is investigated in the process of generating the graded-index profile of plastic optical fibres. Because the diffusion coefficient in high polymers has been shown to depend strongly on dopant concentration, it is allowed in this work to vary with the radial coordinate of the multistep-core fibre. A novel multi-layer model is presented for solving the diffusion equation with the variable diffusion coefficient. It is solved by the finite difference method. The solution determines the dopant diffusion profile in the fibre. It is verified against a solution from the literature and two cases of fibres with diffused profiles. The application is demonstrated on two examples of graded-index plastic optical fibres, one originally with a two-step and the other with four-step core. The results indicate that closer to the core-cladding interface, the computed diffused profile with variable diffusion coefficient D is closer to target profile than the profile obtained with constant D for the same time of thermal process.