A boundary element model provides great flexibility for the simulation of membrane-type micromachined ultrasonic transducers (MUTs) in terms of membrane shape, actuating mechanism, and array layout. Acoustic crosstalk is accounted for through a mutual impedance matrix which captures the primary crosstalk mechanism of dispersive-guided modes generated at the fluid-solid interface. However, finding the solution to the fully-populated boundary element matrix equation using standard techniques requires computation time and memory usage which scales by the cube and by the square of the number of nodes, respectively, limiting simulation to a small number of membranes. We implement a solver with improved speed and efficiency through the application of a multi-level fast multipole algorithm (FMA). By approximating the fields of collections of nodes using multipole expansions of the free-space Green’s function, an FMA solver can enable the simulation of hundreds of thousands of nodes while incurring an approximation error that is controllable. Convergence is drastically improved using a problem-specific block-diagonal preconditioner. We demonstrate the solver’s capabilities by simulating a 32-element 7 MHz 1-D CMUT phased array with 2880 membranes. The array is simulated using 233,280 nodes for a very wide frequency band up to 50 MHz. For a simulation with 15,210 nodes, the FMA solver performed 10-times faster and used 32-times less memory than a standard solver based on LU decomposition We investigate the effects of mesh density and phasing on the predicted array response and find that it is necessary to use about 7 nodes over the width of the membrane to observe convergence of the solution–even below the first membrane resonance frequency–due to the influence of higher-order membrane modes.