  The classical equations of motion for an axion with potential
$V(\phi)=m_a^2f_a^2 [1-\cos (\phi/f_a)]$ possess quasi-stable, localized,
oscillating solutions, which we refer to as "axion stars". We study, for the
first time, collapse of axion stars numerically using the full non-linear
Einstein equations of general relativity and the full non-perturbative cosine
potential. We map regions on an "axion star stability diagram", parameterized
by the initial ADM mass, $M_{\rm ADM}$, and axion decay constant, $f_a$. We
identify three regions of the parameter space: i) long-lived oscillating axion
star solutions, with a base frequency, $m_a$, modulated by self-interactions,
ii) collapse to a BH and iii) complete dispersal due to gravitational cooling
and interactions. We locate the boundaries of these three regions and an
approximate "triple point" $(M_{\rm TP},f_{\rm TP})\sim (2.4 M_{pl}^2/m_a,0.3
M_{pl})$. For $f_a$ below the triple point BH formation proceeds during winding
(in the complex $U(1)$ picture) of the axion field near the dispersal phase.
This could prevent astrophysical BH formation from axion stars with $f_a\ll
M_{pl}$. For larger $f_a\gtrsim f_{\rm TP}$, BH formation occurs through the
stable branch and we estimate the mass ratio of the BH to the stable state at
the phase boundary to be $\mathcal{O}(1)$ within numerical uncertainty. We
discuss the observational relevance of our findings for axion stars as BH
seeds, which are supermassive in the case of ultralight axions. For the QCD
axion, the typical BH mass formed from axion star collapse is $M_{\rm BH}\sim
3.4 (f_a/0.6 M_{pl})^{1.2} M_\odot$.
