The boundary element method (BEM) proves to be a powerful alternative numerical method for solving non-linear diffusion problems. In the BEM the mesh is generated only on the surface of a 3D body or on the contour if the domain is 2D. In order to impose body forces to the problem a technique called as Dual Reciprocity Method (DRM) was introduced. The DRM allows to transform the domain integrals into the boundary equivalent integrals by expanding the inhomogeneous terms into a set of
global approximating basis functions. In this work a numerical routine based on BEM is implemented for solving a non-linear problem and the non-homogeneous terms were dealt by using the DRM. The routine is validated by using a benchmark of a 2D cube model submitted to a thermal load increment
