Understanding the entire history of the ionization state of the intergalactic medium (IGM) is at the frontier of astrophysics and cosmology. A promising method to achieve this is by extracting the damping wing signal from the neutral IGM. As hundreds of redshift z>6 quasars are observed, we anticipate determining the detailed time evolution of the ionization fraction with unprecedented fidelity. However, traditional approaches to parameter inference are not sufficiently accurate.We assess the performance of a simulation-based inference (SBI) method to infer the neutral fraction of the universe from quasar spectra. The SBI adeptly exploit the shape information of the damping wing, enabling precise estimations of the neutral fraction \xHIv and the wing position w_p. Importantly, the SBI framework successfully breaks the degeneracy between these two parameters, offering unbiased estimates of both. This makes the SBI superior to the traditional method using a pseudo-likelihood function. We anticipate that SBI will be essential to determine robustly the ionization history of the Universe through joint inference from the hundreds of high-(z) spectra we will observe.