<p>Non-invasive electromagnetic brain stimulation uses electrodes and/or coils to modulate brain activity via the induced E-fields. E-field dosimetry solvers have improved non-invasive electromagnetic brain stimulation protocol and our understanding of neuroscience. However, E-field dosimetry techniques are incomplete in that the contributions of non-linear neuron activity are left unaccounted for. To better understand the neurological effects of non-invasive electromagnetic stimulation, we introduce an integral equation formulation for modeling the non-linear behavior of neurons due to an incident E-field generated by electrode and coil sources. We formulate the new integral equation using a boundary element approach. We compare the boundary element solver accuracy with an established finite element solver and multi-scale cable equation approaches. Unlike previous approaches, this new boundary integral formulation avoids multi-scaling challenges from meshing while retaining the accuracy and the robust spatial support of integral equation-based methods. The memory savings from switching to surface meshes makes simulations with more complex morphologies computationally tractable. Additionally, we examine the ability of neurons to couple to one another via the local extracellular fields. Examples of simulations with both transcranial electric and magnetic stimulation results for simple geometries are used to illustrate the capabilities of a boundary integral approach. This boundary integral method will aid the development of better neurological understanding, delineate the mechanisms by which electromagnetic stimulation engenders neuronal activity, and aid in modeling local E-field coupling.</p>