This paper develops a method to carry out the large-$N$ asymptotic analysis of a class of $N$-dimensional integrals arising in the context of the so-called quantum separation of variables method. We push further ideas developed in the context of random matrices of size $N$, but in the present problem, two scales $1/N^α$ and $1/N$ naturally occur. In our case, the equilibrium measure is $N^α$-dependent and characterised by means of the solution to a $2\times 2$ Riemann--Hilbert problem, whose large-$N$ behavior is analysed in detail. Combining these results with techniques of concentration of measures and an asymptotic analysis of the Schwinger-Dyson equations at the distributional level, we obtain the large-$N$ behavior of the free energy explicitly up to $o(1)$. The use of distributional Schwinger-Dyson is a novelty that allows us treating sufficiently differentiable interactions and the mixing of scales $1/N^α$ and $1/N$, thus waiving the analyticity assumptions often used in random matrix theory.