  We are interested in the asymptotic behavior of orthogonal polynomials of the
generalized Jacobi type as their degree $n$ goes to $\infty$. These are defined
on the interval $[-1,1]$ with weight function
$w(x)=(1-x)^{\alpha}(1+x)^{\beta}h(x)$, $\alpha,\beta>-1$ and $h(x)$ a real,
analytic and strictly positive function on $[-1,1]$. This information is
available in the work of Kuijlaars, McLaughlin, Van Assche and Vanlessen, where
the authors use the Riemann--Hilbert formulation and the Deift--Zhou non-linear
steepest descent method. We show that computing higher-order terms can be
simplified, leading to their efficient construction. The resulting asymptotic
expansions in every region of the complex plane are implemented both
symbolically and numerically, and the code is made publicly available. The main
advantage of these expansions is that they lead to increasing accuracy for
increasing degree of the polynomials, at a computational cost that is actually
independent of the degree. In contrast, the typical use of the recurrence
relation for orthogonal polynomials in computations leads to a cost that is at
least linear in the degree. Furthermore, the expansions may be used to compute
Gaussian quadrature rules in $\mathcal{O}(n)$ operations, rather than
$\mathcal{O}(n^2)$ based on the recurrence relation.
