We present three-dimensional calculations of spherically symmetric Bondi accretion onto a stationary supermassive black hole (SMBH) of mass $10^{8}$ $M_{odot}$ within a radial range of $0.02-10$ pc, using a modified version of the smoothed particle hydrodynamics (SPH) gad sp code, which ensures approximate first-order consistency (i.e., second-order accuracy) for the particle approximation. First-order consistency is restored by allowing the number of neighbours, $n_{rm neigh}$, and the smoothing length, $h$, to vary with the total number of particles, $N$, such that the asymptotic limits $n_{rm neigh}toinfty$ and $hto 0$ hold as $Ntoinfty$. The ability of the method to reproduce the isothermal ($gamma =1$) and adiabatic ($gamma =5/3$) Bondi accretion is investigated with increased spatial resolution. In particular, for the isothermal models the numerical radial profiles closely match the Bondi solution, except near the accretor, where the density and radial velocity are slightly underestimated. However, as $n_{rm neigh}$ is increased and $h$ is decreased, the calculations approach first-order consistency and the deviations from the Bondi solution decrease. The density and radial velocity profiles for the adiabatic models are qualitatively similar to those for the isothermal Bondi accretion. Steady-state Bondi accretion is reproduced by the highly resolved consistent models with a percent relative error of $lesssim 1$% for $gamma =1$ and $sim 9$% for $gamma =5/3$, with the adiabatic accretion taking longer than the isothermal case to reach steady flow. The performance of the method is assessed by comparing the results with those obtained using the standard Gadget and the Gizmo codes.