For non-linear problems, the solution of the heat equation in terms of the Kirchhoff transformation, θ(T), is very limited. This restriction is due to the practical disadvantage of the inverse temperature shift from the Kirchhoff transform T(θ). In order to get around the difficulties associated with the representation of θ(T) and its inverse T(θ) for solids with strongly non-linear conductivities, a strategy based on a reverse engineering method is considered. It consists in identifying the number of knots and their respective locations on the curve T(θ) at the most efficient computational cost. In order to obtain the location of the knots, the curve is fitted by B-spline functions and the data is partitioned by an application of the bisectional method with a predetermined error. These knots are further optimized using the non-linear least squares method. The proposed approach can be combined with a numerical method such as the FEM, BEM, and FEV to provide the non-linear solution of the heat equation in terms of θ. However, in this work we have limited ourselves to the FEM. The validation of the proposed approach is achieved through several cases ranging from constant to strongly non-linear thermal conductivities with or without convection. As an application, the 3D finite element method is applied to determine the non-linear temperature distribution in a copper block with an imposed temperature.