The boundary element method (BEM) is developed from the standpoint of software design. The Fortran language is used to produce a structured library for solving Laplace's equation in various domain topologies and dimensions with generalised boundary conditions. Subroutines that compute the discrete Laplace operators, which are the core components for populating the matrices in the BEM, are developed. The main subroutines for solving Laplace's equation in 2D, 3D and axisymmetric cases for open and closed boundaries are introduced. The methods are demonstrated on test problems.