A novel numerical method, based on the combined Taylor series method with a technique for rapidly selecting suitable step sizes, is given for solving hva-point boundary value problems, The differential equation for an immobilized enzyme reaction is considered as a model system, Comparisons between numerical and exact solutions for the first-order reaction show that the proposed method promises superhigh-order accuracy that is almost the same as machine accuracy over wide ranges of relevant parameters. The final value of the increment or decrement used to obtain a new estimate in the Newton-Raphson iterative scheme is found to be approximately the same as the relative error of the calculated value, which means that this final value is useful to predict the accuracy of a numerical solution of the differential equation that does not provide an analytical solution. Moreover, it is shown that the selected stepsizes are suitable values that cause little loss-of-significance errors.