An efficient, accurate, and flexible numerical method is proposed for the solution of the swimming problem of one or more autophoretic particles in the purely diffusive limit. The method relies on successive boundary element solutions of the Laplacian and the Stokes flow equations using regularised Green’s functions for swift, simple implementations, an extension of the well-known method of “regularised stokeslets” for Stokes flow problems. The boundary element method is particularly suitable for phoretic problems, since no quantities in the domain bulk are required to compute the swimming velocity. For time-dependent problems, the method requires no re-meshing and simple boundaries such as a plane wall may be added at no increase to the size of the linear system through the method of images. The method is validated against two classical examples for which an analytical or semi-analytical solution is known, a two-sphere system and a Janus particle, and provides a rigorous computational pipeline to address further problems with complex geometry and multiple bodies.