In an earlier work, the neuronal current was expressed via the Helmholtz decomposition in terms of its irrotational part characterised by a scalar function and its solenoidal part characterised by a vectorial function. Furthermore, it was shown that EEG data is affected only by the irrotational part of the current, whereas MEG data is affected by two scalar functions, namely the irrotational component and the radial part of the soleoidal vectorial function. Here, we focus on the numerical implementation of this approach on the three-layer ellipsoidal model. The parameterization of the unknown functions in terms of ellipsoidal harmonics implicitly regularizes the highly ill-posed associated inverse problems. However, despite the above parametrization of these two unknown functions in terms of ellipsoidal harmonics, the inversion matrices are highly ill-conditioned for both EEG and MEG. In order to bypass this problem, we propose an alternative approach to the inversion problem. This involves revisiting the general inversion formulas presented earlier by one of the authors and expressing them as surface integrals. By choosing a suitable parametrization for the relevant unknown functions, these surface integrals can be evaluated using a method for numerical quadrature over smooth, closed surfaces. The method uses local radial basis function interpolation for generating quadrature weights for any given node set. This gives rise to a stable linear system of equations suitable for inversion and reconstruction purposes. We illustrate the effectiveness of our approach by presenting simple reconstructions for both EEG and MEG in a setting where data has a 20 dB SNR.