Standard confidence intervals employed in applied statistical analysis are usually based on asymptotic approximations. Such approximations can be considerably inaccurate in small and moderate sized samples. We derive accurate confidence intervals based on higher-order approximate quantiles of the score function. The coverage approximation error is $O(n^{-3/2})$ while the approximation error of confidence intervals based on the asymptotic normality of MLEs is $O(n^{-1/2})$. Monte Carlo simulations confirm the theoretical findings. An implementation for regression models and real data applications are provided.