quoted from:https://csdms.colorado.edu/wiki/Model:GFlex

The model gFlex implements multiple (user-selectable) solution methods to solve for flexure and isostasy due to surface loading in both one (line loads) and two (point loads) dimensions. It works for elastic lithospheric plates of both constant and spatially variable elastic thickness and allows the user to select the solution method.

An analytical approach to solving the flexure equations is computed by the superimposition of analytical solutions for flexure with a constant elastic thickness. Current implementations perform this superposition in the spatial domain. This works both on uniform grids and arbitrary meshes.

For the numerical implementation, gFlex computes a finite difference solution to the flexure equations for a lithospheric plate of nonuniform (or uniform, if one so desires) elastic thickness via a thin plate assumption. It uses the UMFPACK direct solvers to compute the solutions via a lower-upper decomposition of a coefficient matrix. The coefficient matrix in the 1D case is a pentadiagonal sparse matrix that is trivial to generate. In the 2D case, 2D grid is reordered into a 1D vector. UMFPACK solution routines are then able to copmute solutions to flexure in around a second or less, though sometimes up to a minute for very large grids. Iterative solution methods may also be used, but are not tested.