Abstract We extend the ultraspherical spectral method to solving nonlinear ordinary differential equation (ODE) boundary value problems. Naive ultraspherical Newton implementations usually form dense linear systems explicitly and solve these systems exactly by direct methods, thus suffering from the bottlenecks in both computational complexity and storage demands. Instead, we propose to use the inexact Newton–GMRES framework for which a cheap but effective preconditioner can be constructed and a fast Jacobian-vector multiplication can be effected, thanks to the structured operators of the ultraspherical spectral method. The proposed inexact Newton–GMRES–ultraspherical framework outperforms the naive implementations in both speed and storage, particularly for large-scale problems or problems whose linearization has solution-dependent variable coefficients in higher-order terms. Additional acceleration can be gained when the method is implemented with mixed precision arithmetic.