  We investigate the $\theta$-dependence of 2-dimensional $CP^{N-1}$ models in
the large-$N$ limit by lattice simulations. Thanks to a recent algorithm
proposed by M. Hasenbusch to improve the critical slowing down of topological
modes, combined with simulations at imaginary values of $\theta$, we manage to
determine the vacuum energy density up the sixth order in $\theta$ and up to $N
= 51$. Our results support analytic predictions, which are known up to the
next-to-leading term in $1/N$ for the quadratic term in $\theta$ (topological
susceptibility), and up to the leading term for the quartic coefficient $b_2$.
Moreover, we give a numerical estimate of further terms in the $1/N$ expansion
for both quantities, pointing out that the $1/N$ convergence for the
$\theta$-dependence of this class of models is particularly slow.
