The reconstruction of the neuronal current inside the human brain from magnetic flux density (MEG) and electric potential (EEG) measurements is an important tool for understanding the functioning of the brain and for diagnosing brain diseases, such as epilepsy. One partly unanswered question, which is extensively discussed in the literature, is about the non-uniqueness of the related inverse problems. We investigate this question in the context of the multiple-shell model, which assumes nested spherical geometries. We derive novel integral equations describing the inverse problems, which require less a-priori assumptions on the current than former approaches and map the entire vector-valued current onto the data instead of certain scalar functions. A novel decomposition of the current based on an orthonormal basis system is presented, which yields singular value decompositions with exponentially fast decreasing singular values of the integral operators. Therewith, we complete the existing non-uniqueness considerations, which includes a characterization of the measurable radial part of the neuronal current: only the harmonic solenoidal part of the current can be measured via the MEG and EEG devices. For the numerical solution of these severely ill-posed problems, regularization methods are required. We use the regularized Ritz method, scalar splines, novel vector reproducing kernel based splines, and the regularized (orthogonal) functional matching pursuit (ROFMP) algorithm, which was developed by the Geomathematics Group Siegen within the last years. We improve convergence results of the RFMP and introduce novel Sobolev norms as penalty terms. The good and stable numerical results of the vector spline method are outperformed by the ROFMP in non-noisy and noisy synthetic test cases. Finally, we apply the ROFMP to real measurement data.