An analysis of the eigenmodes associated with a dielectric sphere embedded within one or more spherical layers of varying permittivity is presented. A matrix formulation is applied to the general multilayer problem, and the solution is validated by comparison to the complex transcendental equations derived for the homogeneous dielectric spherical resonator and a single-layer spherical dielectric resonator. Although only lossless cases are considered for simplicity, the formulation presented is applicable to resonators (with loss or gain) via the replacement of real permittivities with complex permittivities. The matrix formulation presented for the multilayer dielectric sphere follows a previously published development with the inclusion of some significant corrections. A Wronskian filter effective index approximation is developed for the estimation of Bessel functions of large order (expanding the functionality of the models to evaluation of whispering gallery modes), which is easily integrated into the presented models. Numerical simulations are generated utilizing MATLAB to solve for the complex roots of the characteristic equation, amplitude reflection coefficient, surface impedance of the dielectric layers, Q-factor, and field patterns for both TE and TM modes. Results are presented for single-layer and multilayer spherical geometries and compared to previously published results.