A theoretical analysis was conducted on a multilayer internal diffusion and surface evaporation (IDSE) coextrusion process for preparing gradient-index (GI) polymer optical fibers (POFs). The predicted refractive index distribution (RID) was in good agreement with the experimental data reported in the literature. The effects of the essential parameters, including the mass transfer coefficient of monomer across solid-gas interface, k, the diffusivities of monomers, and the radius ratio of each layer on RID were investigated. We show that if k is either very small or very large, then RID deviates significantly from a parabolic curve. The effect of the diffusivity of BzMA monomer on RID is more significant than that of MMA monomer. The results of numerical simulation reveal that both the radius ratio and the composition of each layer have significant effect on RID, and the greater the number of layers, the closer the RID is to a parabolic curve.