<p>In this study, we numerically solved the diffusion equation with a concentration-dependent diffusion coefficient. When intrinsic diffusivity depends on the power of concentration, the solution to the diffusion equation was analytically derived in previous studies by Pelton and Etsell (Acta Metall 20: 1269–1274, 1972) and Blanc (J Appl Phys 45: 1948–1950, 1974). We confirmed that analytical methods for solving the diffusion equation are limited to intrinsic diffusion cases where diffusion processes are uncoupled, and we extended the method to thermodynamically coupled interdiffusion cases. The interdiffusion solutions were obtained numerically using a two-step iterative procedure, consisting of an inner step that determines the initial slope of the solution via the shooting method and an outer step that verifies the satisfaction of mass conservation. To validate our extended methodology, solution plots and penetration distances were compared with those reported in previous studies. We analyzed the characteristics of solutions to the diffusion equation where the diffusion coefficient is given as a chemical diffusivity. In particular, we examined the second derivative of the concentration profile under varying the diffusivity. The results show that increased curvature in the diffusion coefficient leads to increased curvature in the concentration profile. This improved our understanding of the nature of diffusion equation solutions in multicomponent systems.</p>