Recent advances in neuronal current imaging using magnetic resonance imaging and in invasive measurement of neuronal magnetic fields have given a need for methods to compute the magnetic field inside a volume conductor due to source currents that are within the conductor. In this work, we derive, verify, and demonstrate an analytical expression for the magnetic field inside an anisotropic multilayer spherically symmetric conductor due to an internal current dipole. We casted an existing solution for electric field to vector spherical harmonic (VSH) form. Next, we wrote an ansatz for the magnetic field using toroidal-poloidal decomposition that uses the same VSHs. Using properties of toroidal and poloidal components and VSHs and applying magnetic scalar potential, we then formulated a series expression for the magnetic field. The convergence of the solution was accelerated by formulating the solution using an addition-subtraction method. We verified the resulting formula against boundary-element method. The verification showed that the formulas and implementation are correct; 99th percentiles of amplitude and angle differences between the solutions were below 0.5% and 0.5 , respectively. As expected, the addition-subtraction model converged faster than the unaccelerated model; close to the source, 250 terms gave relative error below 1%, and the number of needed terms drops fast, as the distance to the source increases. Depending on model conductivities and source position, field patterns inside a layered sphere may differ considerably from those in a homogeneous sphere. In addition to being a practical modeling tool, the derived solution can be used to verify numerical methods, especially finite-element method, inside layered anisotropic conductors. V