Computational geodynamics is a key tool in Earth Sciences, enabling the simulation of different geodynamic processes at any time scale in Nature or those which cannot be adequately reproduced by analogue modelling in laboratorial conditions. Specifically, the coupling of surface processes of erosion and sedimentation with the internal dynamics of the lithosphere is a complex problem that involves the solution of a set of differential equations that are only adequately solved by numerical models. During the last three decades, different numerical models were developed to explain the importance of the coupling between the surface and internal dynamics of the Earth, both in active margins and in stable tectonic domains, showing how the coupling leads to counter-intuitive results not observed when each process is analyzed separately. One example of this complex coupling is the feedback mechanism between erosion of the landscape and the regional isostatic response of the lithosphere. In this work, we present simple isostatic and flexural elements that highlight the importance of surface processes on the stress and strain pattern in lithospheric plates. Initially, we present a review on the development of computational geodynamics at University of São Paulo. This review is followed by an analysis of the density structure of the Earth and how the high-density contrast at the Earth’s surface creates a major impact on the isostatic equilibrium of the lithosphere when variations on topographic loads are taken into account. Additionally, we show that the wavelength of the denudation of the Earth’s surface due to fluvial dynamics corresponds to the characteristic length scale for flexural bending of the lithosphere, maximizing the flexural stresses in the lithosphere. Finally, we present recent works on the coupling between surface and lithospheric processes and future challenges for the development of computational geodynamics, with possible strategies to solve them.