We study Rotating Boson Star initial data for Numerical Relativity as previously considered by Yoshida and Eriguchi, Lai (arXiv:gr-qc/0410040v2), and Grandclement, Some 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 Ruiz, Alcubierre and Nu~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 fields 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.