Considerable progress in the understanding of thin film flow over surfaces has been achieved thanks to lubrication theory which enables the governing Navier-Stokes equations to be reduced to a more tractable form, namely a coupled set of partial differential equations. These are solved numerically since the flows of interest involve substrates containing heterogeneities in the form of wetting patterns and/or topography. An efficient and accurate numerical method is described and used to solve two classes of problem: droplet spreading in the presence of wetting and topographic heterogeneities; gravity-driven flow of continuous thin liquid films down an inclined surface containing well defined topographic features. The method developed, employs a Full Approximation Storage (FAS) multigrid algorithm, is fully implicit and has embedded within it an adaptive time-stepping scheme that enables the same to be optimised in a controlled manner subject to a specific error tolerance. Contact lines are ubiquitous in the context of droplet spreading and the wellknown singularity which occurs there is alleviated by means of a disjoining pressure model. The latter allows prescription of a local equilibrium contact angle and three dimensional numerical simulations reveal how droplets can be forced to either wet or dewet a region containing topography depending on the surface wetting characteristics. The growth of numerical instabilities, in the contact line region, which can lead to the occurrence of non-physical, negative film thicknesses is avoided by using a Positivity Preserving Scheme. A range of two- and three-dimensional problems is explored featuring the gravity-driven flow of a continuous thin liquid film over a non-porous inclined flat surface containing topography. Important new results include: the quantification of the validity range of the lubrication approximation for step-up and step-down topographies; description of the "bow wave" triggered by localised topography and an explanation, in terms of the local flow rate, of the accompanying "downstream surge": an assessment of linear superposition as a means of examining free surface response to topographies. In addition, the potential of local mesh refinement as a means of reducing computational time is highlighted. Finally, more complex liquids composed of a non-volatile resin dissolved in a solvent and allowed to evaporate are considered. An evaporation model based on the wellmixed approximation is utilised. Results show that localised topographies produce defects in dried continuous films which persist far downstream of the topography, while with respect to droplet motion, solvent evaporation is found to be responsible for contact line pinning and thus a reduction in spreading.