Functionally gradient material (FGM) developed for heat-shielding structure is often used in the very high temperature environment. Therefore, the material property parameters are not only functions of spatial coordinates but also ones of temperature. The former leads to partial differential equations with variable coefficients, the latter to nonlinear governing equations. It is usually very difficult to obtain the analytical solution to such thermal conduction problems of FGMs. If the finite element method is adopted, it is very inconvenient because material parameter values must be imputed for each element. Hence, a semi-analytic numerical method, i.e., method of lines (MOLs) is introduced. The thermal conductivity functions do not need to be discretized and remain original form in ordinary differential equations. As a numerical example, the nonlinear steady temperature fields are computed for a rectangular non-homogeneous region with the first, the second and the third kinds of boundary conditions, where three kinds of functions, i.e. power, exponential and logarithmic ones are adopted for the thermal conductivity. Results display the important influence of non-linearity on temperature fields. Moreover, the results given here provide the better basis for thermal stress analysis of non-homogenous and non-linear materials.