In this contribution, we modeled the thermal field of the colombian lithosphere by solving the 2D heat equation along three profiles parallel (1) and perpendicular (3) to the Eastern Cordillera. For this purpose, the conduction-advection equation was discretized by the conservative finite difference method. This numerical solution was validated by comparing results with analytical solutions of some simple 1D and 2D tectonic scenarios reported in the literature. We use this numerical approach to solve the lithospheric thermal field in the Eastern Cordillera of Colombia. The Geological Survey of Colombia geothermal gradient database and some heat flow estimates in the Eastern Llanos basin were used to constrain the numerical solution. We assigned typical values to parameters such as thermal conductivity and radiogenic heat production, among others, to rocks that were outcropped in this region and inferred their distribution in the crust and upper mantle by global thickness models. The boundary conditions and values of parameters that best fit the surface data suggest the presence of a basement enriched in radioactive elements that would be an important source of heat (1.6 – 3 μWm-3) at the east of the Guaicaramo Thrust Fault System, where active geothermal systems have been documented. The modeled asthenosphere heat flow ranges in 15-25 mWm-2 except on the Eastern Cordillera, north of the Caldas tear, where heat flow needed to raise up the lithosphere-asthenosphere boundary to 70-80 km, requires higher values (30-45 mWm-2). The additional heat flow could be due to advection by mantle flow not considered in the simulation or by viscous heating in the mantle wedge.