  We present an algorithm for constructing numerical solutions to
one--dimensional nonlinear, variable coefficient boundary value problems. This
scheme is based upon applying the Homotopy Analysis Method (HAM) to decompose a
nonlinear differential equation into a series of linear differential equations
that can be solved using a sparse, spectrally accurate Gegenbauer
discretisation. Uniquely for nonlinear methods, our scheme involves
constructing a single, sparse matrix operator that is repeatedly solved in
order to solve the full nonlinear problem. As such, the resulting scheme scales
quasi-linearly with respect to the grid resolution. We demonstrate the
accuracy, and computational scaling of this method by examining a fourth-order
nonlinear variable coefficient boundary value problem by comparing the scheme
to Newton-Iteration and the Spectral Homotopy Analysis Method, which is the
most commonly used implementation of the HAM.
