Let $P_{n}(z)=\gamma_{n}z^{n}+\cdots$ be a sequence of orthonormal polynomials with respect to some weight $w$ on $\mathbb{R}$. Given an $n\geq0$, consider the following Riemann-Hilbert problem for the $2\times2$ matrix-valued function $Y_{n}(z)$ of the complex variable $z$,
1) $z\to Y_{n}(z)$ is analytic on $\mathbb{C}\setminus\mathbb{R}$,
2) $Y_{n,+}(z)=Y_{n,-}(z)\begin{pmatrix} 1 & w\\0 & 1\end{pmatrix}$, where $Y_{n,\pm}$ denotes the upper and lower limits of $Y_{n}$ on $\mathbb{R}$,
3) $Y_{n}(z)=(I+\mathcal{O}(1/z))\begin{pmatrix} z^{n} & 0 \\ 0 & z^{-n}\end{pmatrix}$ as $z\to\infty$.
The Riemann-Hilbert problem has a unique solution $Y_{n}$, which is given by
$$Y_{n}(z)=\begin{pmatrix} P_{n}(z) & {\displaystyle\dfrac{1}{2i\pi}\int\dfrac{P_{n}(x)w(x)}{x-z}dx}\\[10pt]
c_{n-1}P_{n-1}(z) & {\displaystyle\dfrac{c_{n-1}}{2i\pi}\int\dfrac{P_{n-1}(x)w(x)}{x-z}dx}
\end{pmatrix},
\qquad c_{n-1}=-2i\pi\gamma_{n-1}^{2},$$
(the jump of $Y_n$ on the real line comes from the singular integrals on the second row).
Performing a series of transformations on the initial Riemann-Hilbert problem leads to a Riemann-Hilbert problem normalized at infinity, and whose jump matrix is uniformly close to the identity matrix as $n$ tends to infinity. Using the asymptotic expansion in powers of $1/n$ of this jump matrix, one may derive an asymptotic expansion for $Y_{n}$ and thus for $P_{n}$. This method is described in
the book ``Orthogonal polynomials and random matrices: a Riemann-Hilbert approach'' by P. Deift.
As is well-known, the polynomials $P_{n}$ satisfy a three-term recurrence relation
$$
b_{n-1}P_{n-1}(z)+(a_{n}-z)P_{n}(z)+b_{n}P_{n+1}(z)=0,
$$
and asymptotic expansions of the coefficients $a_{n}$ and $b_{n-1}$ can also be derived by the above method because
$$
a_{n}=(Y_{n,1})_{11}-(Y_{n+1,1})_{11},\qquad b_{n-1}^{2}=(Y_{n,1})_{12}(Y_{n,1})_{21},
$$
where
$$
Y_{n}(z)=(I+z^{-1}Y_{n,1}+\mathcal{O}(z^{-2}))\begin{pmatrix}
z^{n} & 0\\ 0 & z^{-n}
\end{pmatrix},
$$
see Chapter 3 of Deift's book. Finally, the coefficients in the recurrence relation also appear in the expansion as a continued fraction of
$$
\frac{Q_{n}(z)}{P_{n}(z)}=\frac{1}{z-a_{0}-{\displaystyle\frac{b_{0}^{2}}{z-a_{1}-\cdots}}},
$$
where $Q_{n}$ is the associated polynomial of the second kind. For an application of the method to polynomials orthogonal with respect to the modified Jacobi weight $(1-x)^{\alpha}(1+x)^{\beta}h(x)$ on $[-1,1]$, $h$ a positive real analytic function, see Arxiv, where, for instance, asymptotic expansions for the recurrence coefficients are explicitly computed (Theorem 1.10).
