We present a detailed phenomenological study of forward hadron neutral pion production in deep inelastic scattering, with both the direct and the resolved contributions calculated to NLO accuracy. A comparison of the theoretical predictions for the various distributions with the H1 data and a study of stability of the QCD predictions under changes of scales is the focus of this study. We obtain a very good overall description of the recent H1 data with the choice of scale Q^2 + E^2_{perp}, in contrast to the (Q^2 + E^2_{perp})/2 required earlier when the resolved contribution was included only at LO accuracy. We find a more modest variation of the predictions, as the scale is changed from (Q^2 + E^2_{perp})/2 to 2(Q^2 + E^2_{perp}), as compared to the case where the resolved contribution was included only at LO accuracy. This variation is of the order of the rather large experimental errors. Unfortunately, this fact prevents us from concluding that perturbation theory gives an unambiguous prediction for forward particle production in deep inelastic scattering. However, the overall success of perturbative QCD in explaining the small x_{Bj} data means that perhaps a full resummation of the BFKL ladder is not called for. We notice the need for rather large resolved contributions to explain the data at low x_{Bj} even at somewhat larger Q^2 values.