Thin plate flexure theory provides an accurate model for the response of the lithosphere to vertical loads on horizontal length scales ranging from tens to hundreds of kilometres. Examples include flexure at seamounts, fracture zones, sedimentary basins and subduction zones. When applying this theory to real world situations, most studies assume a locally uniform plate thickness to enable simple Fourier transform solutions. However, in cases where the amplitude of the flexure is prominent, such as subduction zones, or there are rapid variations in sea floorage, such as fracture zones, these models are in adequate. Here we present a computationally efficient algorithm for solving the thin plate flexure equation for non-uniform plate thickness and arbitrary vertical load. The iteratives cheme takes advantage of the 2-D fast Fourier transform to perform calculations in both the spatial and spectral domains, resulting in an accurate and computationally efficient solution. We illustrate the accuracy of the method through comparisons with known analytic solutions. Finally, we present results from three simple models demonstrating the differences in trench outer rise flexure when 2-D variations in plate rigidity and applied bending moment are taken into account. Although we focus our analysis on ocean trench flexure, theme thodisapplicable too ther 2-D flexure problems having spatialr igidity variations such as sea mount loading of a thermally eroded lithosphere or flexure across the continental–oceanic crust boundary.