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_all M_{pl}$. For larger $f_agtrsim 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$.