  We present a novel numerical solver for the systems of coupled non-linear
elliptical differential equations. The solver partitions the computational
domain into a set of rectangular pseudo-spectral collocation subdomains and is
especially well-suited for working with stiff solutions such as almost
shell-like solitonic boson stars. The method can be used in any number of
dimensions although is practical for one- to three-dimensional problems. We
apply the method to rotating and spherically symmetric solitonic boson stars
and demonstrate that it displays exponential convergence. In the spherical
symmetric case we explore families of almost shell-like solitonic boson stars
and get the results that conform with the well-known analytic approximation.
