  We study Rotating Boson Star initial data for Numerical Relativity as
previously considered by Yoshida and Eriguchi, Lai (arXiv:gr-qc/0410040v2), and
Grandclement, Som\'e and Gourgoulhon (arXiv:1405.4837v3). We use a 3 + 1
decomposition as presented by Gourgoulhon (arXiv:1003.5015v2) and Alcubierre,
adapted to an axisymmetric quasi-isotropic spacetime with added regularization
at the axis following work by Ru\'iz, Alcubierre and N\'u\~nez
(arXiv:0706.0923v2) and Torres. The Einstein-Klein-Gordon equations result in a
system of six-coupled, elliptic, nonlinear equations with an added unknown for
the scalar field's frequency $\omega$. Utilizing a Cartesian two-dimensional
grid, finite differences, Global Newton Methods adapted from Deuflhard, the
sparse direct linear solver PARDISO, and properly constraining all variables
generates data sets for rotation azimuthal integers $l \in [0, 6]$. Our
numerical implementation, published in GitHub, is shown to correctly converge
both with respect to the resolution size and boundary extension (fourth-order
and third-order, respectively). Thus, global parameters such as the Komar
masses and angular momenta can be precisely calculated to characterize these
spacetimes. Furthermore, analyzing the full family at fixed rotation integer
produces maximum masses and minimum frequencies. These coincide with previous
results in literature for $l \in [0,2]$ and are new for $l > 2$. In particular,
the study of high-amplitude and localized scalar fields in axial symmetry is
revealed to be only possible by adding the sixth regularization variable.
