The carbon diffusion in steel, where the carbon diffusivity varies with the carbon content, was solved with the integral methods under the third boundary condition. The variation of carbon diffusivity in steel with the carbon content was described with two different functions, linear dependence and exponential dependence. The integral approximation for both cases was improved with the numerical computation to more accurately predict the carbon profiles. The integral solution is more accurate than the formulation based on the assumption of a constant diffusivity or those based on the assumption of a constant diffusivity and/or constant carbon content at part surface. It is also more easily used in practice than the numerical method to describe the carburising process and predict the carbon content at steel surface and carbon profiles in treated layer.