We study the Bose-Einstein condensate trapped in a three-dimensional spherically symmetrical potential. Exact solutions to the stationary Gross-Pitaevskii equation are obtained for properly modulated radial nonlinearity. The solutions contain vortices with different winding numbers and exhibit the shell-soliton feature in the radial distributions. [8,9] or in space [10,14]. It has been used to generate bright solitons [11,12] and induce collapse [13]. In addition, bright solitons [15] and periodic wave solutions [16] were obtained in spinor BECs. However, only one-dimensional matter-wave solitons have been created in experiments by using cigar-shaped traps to confine the condensate [17]. The realization of higher-dimensional matter-wave solitons is still a challengeable task because they are usually unstable for the nonlinear Schrödinger (NLS) equation with constant or uniform couplings due to the weak and strong collapses of the BECs [7]. Modulation of atomic scattering length by the Feshbach resonance is expected to dynamically stabilize higher dimensional bright solitons [18].On the other hand, a multi-quantized vortex is usually unstable when the BECs are trapped in harmonic potential [19][20][21][22][23]. It will split into several singly-quantized vortices. However, in the presence of a plug potential [24], an anharmonic trapping potential [25,26], or an immiscible two-species BEC [4,5], the multi-quantized vortex can be effectively stabilized. For example, when the confining potential is steeper than the harmonic potential, multi-quantized vortices are energetically favorable.In this paper, we study 3D solutions to the Gross-Pitaevskii equation (GPE) by means of the similar transformation. Exact solitary vortices in the BEC trapped in an external potential are constructed by properly modulating the radial nonlinearity. The number of vortex-soliton (VS) modes is specified by the radial nodes in the wavefunction and the winding number of the vortices.The scaled stationary GPE for the macroscopic wave function readswhere V (r) is the spherically harmonic oscillator potential. g(r) is the nonlinearity coefficient which is spatially variable. We rewrite the wavefunction in the spherical coordinates aswhere Y lm (θ, ϕ) = P |m| l (cos θ)e imϕ (l = 0, 1, 2, · · · , and m = 0, ±1, ±2, · · · , ±l) is the spherical harmonic function. The radial part Φ(r) obeys the following nonlinear equation,The convergence condition for r → 0 requires that Φ(r) ∼ r l for l = 0 while Φ ′ (r) = 0 for l = 0. In the meantime, Φ(r → ∞) → 0 due to the localization of the BEC trapped in the potential.According to the similar transformation, we define Φ(r) ≡ ρ(r)U [R(r)] so as to obtain two independent equations[18],