A method was developed for efficient boundary element method solution of potential problems that were large, linear and involved optimization of boundary geometry. The method involved forming an optimization region containing the section of boundary to be optimized, and an invariant region incorporating the remainder of the problem. Solution was by Gaussian elimination. Matrix elements arising from boundary discretization and application of the governing equations were manipulated to permit computations relating to the invariant region to be performed only once. The matrix and vector terms for the invariant region were first calculated, the matrix converted to upper-diagonal form and the values stored. For each boundary configuration within the optimization region, additional matrix elements were computed and appended to the pre-computed matrix. The complete set of equations were subsequently solved by conversion of the appended matrix terms to upper-diagonal form, followed by back-substitution. A practical example involved determining the electric field arising due to different electrodes implanted in the human heart. Computation time was found to be an order of magnitude less for the second and subsequent electrode configurations than for the first. The method may find application in boundary element method optimization problems involving complex or large regions.