We present a novel finite difference implementation of a three-dimensional higher-order ice sheet model that performs well both in terms of convergence rate and numerical stability. In order to achieve these benefits the discretisation of the governing force balance equation makes extensive use of information on staggered grid points. Using the same iterative solver, an existing discretisation that operates exclusively on the regular grid serves as a reference. Participation in the ISMIP-HOM benchmark indicates that both discretisations are capable of reproducing the higher-order model inter-comparison results. This allows a direct comparison not only of the resultant velocity fields but also of the solver's convergence behaviour which holds main differences. First and foremost, the new finite difference scheme facilitates convergence by a factor of up to 7 and 2.6 in average. In addition to this decrease in computational costs, the precision for the resultant velocity field can be chosen higher in the novel finite difference implementation. For high precisions, the old discretisation experiences difficulties to converge due to large variation in the velocity fields of consecutive Picard iterations. Finally, changing discretisation prevents build-up of local field irregularites that occasionally cause divergence of the solution for the reference discretisation. The improved behaviour makes the new discretisation more reliable for extensive application to real ice geometries. Higher precision and robust numerics are crucial in time dependent applications since numerical oscillations in the velocity field of subsequent time steps are attenuated and divergence of the solution is prevented. Transient applications also benefit from the increased computational efficiency.